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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedSeg_tool.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-tools/tinatool/tlmedical/tlmedSeg_tool.c,v $
 23  * Date    :  $Date: 2008/12/02 22:04:19 $
 24  * Version :  $Revision: 1.10 $
 25  * CVS Id  :  $Id: tlmedSeg_tool.c,v 1.10 2008/12/02 22:04:19 paul Exp $
 26  *
 27  * Notes :
 28  *
 29  *  Gio Buonaccorsi commented out prob_all_plot_proc as it is required only for testing.
 30  *  19 Feb 2003.
 31  *
 32  *********
 33 */
 34 
 35 #include "tlmedSeg_tool.h"
 36 
 37 
 38 #if HAVE_CONFIG_H
 39 #   include <config.h>
 40 #endif
 41 
 42 #include <stdio.h>
 43 #include <math.h>
 44 #include <tina/sys/sysDef.h>
 45 #include <tina/sys/sysPro.h>
 46 #include <tina/math/mathDef.h>
 47 #include <tina/math/mathPro.h>
 48 #include <tina/file/filePro.h>
 49 #include <tina/image/imgDef.h>
 50 #include <tina/image/imgPro.h>
 51 #include <tinatool/draw/drawDef.h>
 52 #include <tinatool/draw/drawPro.h>
 53 #include <tinatool/wdgts/wdgtsPro.h>
 54 #include <tinatool/tlbase/tlbasePro.h>
 55 #include <tinatool/tlmedical/tlmedSeg_hist.h>
 56 
 57 #define MAXPATHLEN 1024
 58 #define GMAX 7.0
 59 #define MAXPLOT 255
 60 
 61 extern FILE *fopen_2();
 62 
 63 
 64 static Mixmodel *model = NULL, *model_2d=NULL;
 65 static void ***pim_array=NULL;
 66 void   *p_em_str_name = NULL;
 67 void   *pfile = NULL;
 68 static char em_str_name[MAXPATHLEN];
 69 static char file_name[MAXPATHLEN];
 70 static int prob_flag = -1;
 71 
 72 static void mD_slope_proc()
 73 {
 74   Imrect *im_slope = NULL;
 75   Sequence *seq = NULL;
 76 
 77   seq = seq_get_current();
 78   if (seq == NULL) return;
 79 
 80   im_slope = (Imrect *)mD_slope_image(seq);
 81   stack_push(im_slope, IMRECT, im_free);
 82   imcalc_draw(imcalc_tv_get());
 83 }
 84 
 85 static void mD_grad_scat_proc()
 86 {
 87    Imrect *im=NULL, *im_slope=NULL, *im_scat=NULL;
 88    Imrect *im_slope2;
 89    Imregion *roi;
 90    List *im_ptr;
 91    Sequence *seq = seq_get_current();
 92 
 93    if (seq == NULL)
 94    {
 95       error("Sequence is missing!", warning);
 96       return;
 97    }
 98 
 99    im_ptr = get_seq_current_el(seq);
100    if(im_ptr->to == NULL)
101    {
102       error("There is no image in sequence",warning);
103       return;
104    }
105    im = im_ptr->to;
106 
107    im_slope = (Imrect *)mD_slope_image(seq);
108    if ((roi = tv_get_im_roi(seq_tv_get())) != NULL)
109    {
110         im_slope2 = im_subim(im_slope, roi);
111         im_free(im_slope);
112         im_slope = im_slope2;
113    }
114    //im_scat = (Imrect *) scat_grad(im,im_slope);
115 
116    im_scat = (Imrect *) mD_scat_grad(im,im_slope);
117    stack_push((void *)im_scat, IMRECT, im_free);
118    imcalc_draw(imcalc_tv_get());
119    im_free(im_slope);
120 }
121 
122 static void md_grad_model_proc()
123 {
124    Mixmodel *m=em_get_model(), *m_sub;
125    Imrect *im=NULL, *im_slope=NULL, *im_model_norm=NULL;
126    List *im_ptr;
127    Sequence *seq = seq_get_current();
128    int im_no;
129 
130    if (seq == NULL)
131    {
132       error("There is no image sequence!",warning);
133       return;
134    }
135    if (m == NULL)
136    {
137        error("Input file is missing!", warning);
138        return;
139    }
140    im_ptr = get_seq_current_el(seq);
141    if(im_ptr->to == NULL)
142    {
143       error("There is no image in sequence!",warning);
144       return;
145    }
146    im = im_ptr->to;
147    im_no = get_seq_current_frame(seq);
148    m_sub = get_submodel(m,im_no,im_no);
149 
150    im_slope = (Imrect *)mD_slope_image(seq);
151    im_model_norm =
152    (Imrect *)mD_grad_model_norm(im,im_slope,m,m_sub);
153    stack_push((void *)im_model_norm, IMRECT, im_free);
154    imcalc_draw(imcalc_tv_get());
155    im_free(im_slope);
156 }
157 
158 void kgrad_calc_proc()
159 {
160    Mixmodel *m=em_get_model();
161    Imrect *im=NULL, *im_slope=NULL, *im_slope2, *im_scat=NULL;
162    Imregion *roi;
163    List *im_ptr;
164    Sequence *seq = seq_get_current();
165    int im_no;
166 
167    if (seq == NULL)
168    {
169       error("There is no image sequence!",warning);
170       return;
171    }
172    if (m == NULL)
173    {
174        error("Input file is missing!", warning);
175        return;
176    }
177 
178    im_ptr = get_seq_current_el(seq);
179    if(im_ptr->to == NULL)
180    {
181       error("There is no image in sequence!",warning);
182       return;
183    }
184    im = im_ptr->to;
185    im_no = get_seq_current_frame(seq);
186 
187    im_slope = (Imrect *)mD_slope_image(seq);
188    if ((roi = tv_get_im_roi(seq_tv_get())) != NULL)
189    {
190         im_slope2 = im_subim(im_slope, roi);
191         im_free(im_slope);
192         im_slope = im_slope2;
193    }
194    im_scat = (Imrect *) mD_scat_grad(im,im_slope);
195    m->gradscale = mD_grad_kcalc(m,im,im_slope,pim_array,NULL);
196 
197    im_free(im_slope);
198    im_free(im_scat);
199 }
200 
201 
202 int *get_prob_flag()
203 {
204    return(&prob_flag);
205 }
206 
207 Mixmodel *em_get_model(void)
208 {
209     return(model);
210 }
211 
212 /* Function to set the mixmodel: acts as a wrapper for calling */
213 /* input_mixmodel from outside this function     PAB 1/1/2005  */
214 
215 void em_set_model(Mixmodel *model_to_set)
216 {
217    int i, j;
218 
219    if (pim_array != NULL && model != 0)
220    {
221       for (i=0; i < model->nmix ; i++)
222          for (j=0; j< model->nmix; j++)
223             im_free(pim_array[i][j]);
224       parray_free(pim_array,0,0,model->nmix,model->nmix);
225    }
226    model = model_to_set;
227 }
228 
229 char *em_get_str(void)
230 {
231    return(em_str_name);
232 }
233 
234 void my_log_proc(void) 
235 {
236    FILE   *stream = NULL;
237    double tmp;
238    int i,j,k;
239 
240    stream = (FILE *)fopen_2("log.txt","a");
241    if (stream == NULL) return;
242 
243    fprintf(stream,"%d %d %f \n",model->nmix,model->ndim,model->background);
244    for (i=0;i<model->nmix;i++)
245    {
246       for (j=0; j<model->ndim; j++)
247       {
248          fprintf(stream," %f",model->vectors[i][j]);
249       }
250      fprintf(stream," %f \n", model->volume[i]);
251    }
252    for (k=0; k< model->nmix; k++)
253    {
254       for (i=0;i<model->ndim;i++)
255       {
256          for (j=0; j<model->ndim; j++)
257          {
258             tmp = model->cov[k][i][j];
259             if (tmp!=0.0) tmp /= sqrt(fabs(tmp));
260             fprintf(stream," %f",tmp);
261          }
262          fprintf(stream, " \n");
263       }
264    }
265    for (i=0;i<model->nmix;i++)
266    {
267       for (j=0; j<model->nmix; j++)
268       {
269          fprintf(stream," %f",model->par_dens[i][j]);
270       }
271       fprintf(stream, " \n");
272    }
273    fclose(stream);
274 }
275 
276 static void ptest_proc1()
277 {
278    void ***imptrs = NULL;
279    int   imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
280    Imrect *im_model=NULL, *im_data=NULL;
281    float *row=NULL, *data=NULL;
282    double **a,**at,a_all=0.0,b;
283    int lx, ux, ly, uy;
284    int width=256, height=256;
285    int i, j, k, loop, loop1, cx,cy;
286    double grey;
287    float xy[2];
288    Sequence *seq = seq_get_current();
289    List *lptr1;
290    int im1_no, im2_no;
291     Imrect *im_1 = NULL, *im_2 = NULL;
292    Imregion *roi=NULL;
293    Imrect *im_seq =  seq_image_get();
294    float high=256, low1=0.0, low2=0.0, min1, max1, scale1, min2, max2,scale2;
295 
296    if (im_seq == NULL)
297    {
298       error("Image is missing!", warning);
299       return;
300    }
301 
302    if (model == NULL)
303    {
304        error("Input file is missing!", warning);
305        return;
306    }
307    lptr1 = get_seq_current_el(seq);
308    if(lptr1->to == NULL || lptr1->next == NULL)
309    {
310       error("There are not two images in sequence", warning);
311       return;
312    }
313    im_1 = lptr1->to;
314    im_2 = lptr1->next->to;
315    im1_no = get_seq_current_frame(seq);
316    im2_no = im1_no+1 ;
317    format(" im1_no = %i im2_no = %i\n",im1_no,im2_no);
318 
319    roi = im_1->region;
320       if (roi == NULL)  return;
321 
322    if(model->ndim >= 2 && im2_no > im1_no)
323       model_2d = get_submodel(model, im1_no,im2_no);
324    else
325    {
326       error("Insufficient images in sequence", warning);
327       return;
328    }
329 
330    a  = darray_alloc(0,0,model->nmix,model->nmix);
331    at = darray_alloc(0,0,model->nmix,model->nmix);
332 
333    imf_minmax(im_1, &min1, &max1);
334    if (max1 == min1)
335       max1 = (float) (min1 + 1.0);
336    scale1 = (float) fabs((high - low1) / (max1 - min1));
337    low1 -= (scale1 * min1);
338 
339    imf_minmax(im_2, &min2, &max2);
340    if (max2 == min2)
341       max2 = (float) (min2 + 1.0);
342    scale2 = (float)fabs((high - low2) / (max2 - min2));
343    low2 -= (scale2 * min2);
344 
345 
346    im_model = im_alloc(height, width, roi, float_v);
347    im_data  = im_alloc(height, width, roi, float_v);
348    lx = roi->lx;
349    ux = roi->ux;
350    ly = roi->ly;
351    uy = roi->uy;
352    format(" lx = %i ux = %i ly = %i uy = %i\n",lx,ux,ly,uy);
353 
354    row = fvector_alloc(lx, ux);
355 
356    for (i = ly; i < uy; ++i)
357    {
358       for (j = lx; j < ux; ++j)
359       {
360          xy[0] = (j  - low1)/scale1;
361          xy[1] = (i  - low2)/scale2;
362          for(loop = 0; loop < model->nmix; loop++)
363          {
364             for(loop1 = 0; loop1 < model->nmix; loop1++)
365             {
366                if (loop == loop1)
367                   a[loop][loop1]=model->normal[loop]*pure_density(model_2d,xy,loop);
368 
369                else
370                {
371                  b = part_fraction(model_2d,xy,loop,loop1);
372                  a[loop][loop1]=model->par_dens[loop][loop1]*
373                                 part_density(model_2d,xy,loop,loop1,b);
374                }
375 
376                a_all += a[loop][loop1];
377                at[loop][loop1]+=1/(scale1*scale2)*a[loop][loop1];
378              }
379           }
380           row[j] = a_all;
381           a_all = 0.0;
382        }
383        im_put_rowf(row,im_model,i, lx, ux);
384     }
385 
386     for(loop = 0; loop < model->nmix; loop++)
387     {
388       for(loop1 = 0; loop1 < model->nmix; loop1++)
389          format(" %f ",at[loop][loop1]);
390       format("\n");
391     }
392     stack_push((void *)im_model, IMRECT, im_free);
393     if (imcalc_tv_get()!=NULL) tv_init(imcalc_tv_get());
394     if (imcal2_tv_get()!=NULL) tv_init(imcal2_tv_get());
395     image_choice_reset();
396     fvector_free(row, lx);
397     darray_free(a,0,0,model->nmix,model->nmix);
398     darray_free(at,0,0,model->nmix,model->nmix);
399 
400     seq_slice_init(seq);
401     imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
402                          &imptrux);
403     data = (float *)fvector_alloc(imptrlz, imptruz);
404 
405     cx = imptrlz + im1_no;
406     cy = cx + 1;
407     for (i = ly; i < uy; ++i)
408     {
409         for (j = lx; j < ux; ++j)
410         {
411            for (k = cx; k <= cy; k++)
412            {
413               row = imptrs[k][i];
414               data[k] = row[j];
415            }
416            grey = 1.0 + im_get_pixf(im_data,(int)(data[cy]*scale2+low2),
417                                             (int)(data[cx]*scale1+low1));
418            im_put_pixf(grey,im_data,(int)(data[cy]*scale2+low2),
419                                     (int)(data[cx]*scale1+low1));               }
420     }
421     stack_push((void *)im_data, IMRECT, im_free);
422     if (imcalc_tv_get()!=NULL) tv_init(imcalc_tv_get());
423     if (imcal2_tv_get()!=NULL) tv_init(imcal2_tv_get());
424     image_choice_reset();
425     fvector_free(data,imptrlz);
426     mixmodel_free(model_2d);
427 }
428 
429 static void ptest_proc()
430 {
431    void ***imptrs = NULL;
432    int   imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
433    Imrect *im_model, *im_data;
434    float *row, *data;
435    double **a,**at,a_all=0.0,b;
436    int lx, ux, ly, uy;
437    int width=256, height=256;
438    int i, j, k, loop, loop1,im1_min,im1_max,im2_min,im2_max, cx,cy;
439    double grey;
440    float xy[2];
441    Sequence *seq = seq_get_current();
442    List *lptr1;
443    int im1_no, im2_no,xmax,ymax,xmin,ymin,max_g, min_g;
444    float scale_const,low=0.0;
445    Imrect *im_1 = NULL, *im_2 = NULL;
446    Imregion roi;
447    Imrect *im_seq =  seq_image_get();
448    Imregion *roi1;
449 
450    if (im_seq == NULL)
451    {
452       error("Image is missing!", warning);
453       return;
454    }
455    roi1 = im_seq->region;
456    if (roi1 == NULL)
457         return;
458    if (model == NULL)
459    {
460        error("Input file is missing!", warning);
461        return;
462    }
463    lptr1 = get_seq_current_el(seq);
464    if(lptr1->to == NULL || lptr1->next == NULL)
465    {
466       error("There are not two images in sequence", warning);
467       return;
468    }
469    im_1 = lptr1->to;
470    im_2 = lptr1->next->to;
471    im1_no = get_seq_current_frame(seq);
472    im2_no = im1_no+1 ;
473    format(" im1_no = %i im2_no = %i\n",im1_no,im2_no);
474 
475    if(model->ndim >= 2 && im2_no > im1_no)
476       model_2d = get_submodel(model, im1_no,im2_no);
477    else
478    {
479       error("Insufficient images in sequence", warning);
480       return;
481    }
482 
483    a  = darray_alloc(0,0,model->nmix,model->nmix);
484    at = darray_alloc(0,0,model->nmix,model->nmix);
485 
486    im1_min = (int)imf_min(im_1, &ymin,&xmin);
487    im1_max = (int)imf_max(im_1, &ymax,&xmax);
488    im2_min=  (int)imf_min(im_2, &ymin,&xmin);
489    im2_max = (int)imf_max(im_2, &ymax,&xmax);
490 
491    format(" im1_min = %i im1_max = %i im2_min = %i im2_max =%i\n",im1_min,im1_max,im2_min,im2_max);
492    if (im1_max> im2_max ) max_g = im1_max;
493    else max_g = im2_max;
494    if (im1_min < im2_min ) min_g = im1_min;
495    else min_g = im2_min;
496    scale_const = fabs(max_g - min_g)/width;
497    if (scale_const < 1) scale_const = 1.0;
498    format(" scale_c = %f max_g = %i min_g = %i\n",scale_const, max_g,min_g);
499    low -= min_g/scale_const;
500 
501    roi.lx = (int) im1_min/scale_const;
502    roi.ux = (int) im1_max/scale_const;
503    roi.ly = (int) im2_min/scale_const;
504    roi.uy = (int) im2_max/scale_const;
505 
506    im_model = im_alloc(height, width, &roi, float_v);
507    im_data  = im_alloc(height, width, &roi, float_v);
508    lx = roi.lx;
509    ux = roi.ux;
510    ly = roi.ly;
511    uy = roi.uy;
512    format(" lx = %i ux = %i ly = %i uy = %i\n",lx,ux,ly,uy);
513 
514    row = fvector_alloc(lx, ux);
515 
516    for (i = ly; i < uy; ++i)
517    {
518       for (j = lx; j < ux; ++j)
519       {
520          xy[0] = scale_const*(j+0.5);
521          xy[1] = scale_const*(i+0.5);
522          for(loop = 0; loop < model->nmix; loop++)
523          {
524             for(loop1 = 0; loop1 < model->nmix; loop1++)
525             {
526                if (loop == loop1)
527                   a[loop][loop1]=model->normal[loop]*pure_density(model_2d,xy,loop);
528                else
529                {
530                   b = part_fraction(model_2d,xy,loop,loop1);
531                   a[loop][loop1]=model->par_dens[loop][loop1]*
532                                  part_density(model_2d,xy,loop,loop1,b);
533                }
534                  //a[loop][loop1] = 0.0;
535                a_all += a[loop][loop1];
536                at[loop][loop1]+=scale_const*scale_const*a[loop][loop1];
537              }
538           }
539           row[j] = a_all;
540           a_all = 0.0;
541        }
542        im_put_rowf(row,im_model,i, lx, ux);
543     }
544     for(loop = 0; loop < model->nmix; loop++)
545     {
546       for(loop1 = 0; loop1 < model->nmix; loop1++)
547          format(" %f ",at[loop][loop1]);
548       format("\n");
549     }
550     stack_push((void *) im_model, IMRECT, im_free);
551     tv_init(imcalc_tv_get());
552     tv_init(imcal2_tv_get());
553     image_choice_reset();
554     fvector_free(row, lx);
555 
556     darray_free(a,0,0,model->nmix,model->nmix);
557     darray_free(at,0,0,model->nmix,model->nmix);
558 
559     lx = roi1->lx;
560     ux = roi1->ux;
561     ly = roi1->ly;
562     uy = roi1->uy;
563 
564     format(" lx = %i ux = %i ly = %i uy = %i\n",lx,ux,ly,uy);
565     seq_slice_init(seq);
566     imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
567                          &imptrux);
568     data = (float *)fvector_alloc(imptrlz, imptruz);
569 
570     cx = imptrlz + im1_no;
571     cy = cx + 1;
572     for (i = ly; i < uy; ++i)
573     {
574         for (j = lx; j < ux; ++j)
575         {
576            for (k = cx; k <= cy; k++)
577            {
578               row = imptrs[k][i];
579               data[k] = row[j];
580            }
581            grey = 1.0 + im_get_pixf(im_data,(int)(data[cy]/scale_const),
582                                             (int)(data[cx]/scale_const));
583            im_put_pixf(grey,im_data,(int)(data[cy]/scale_const),
584                               (int)(data[cx]/scale_const));         }
585     }
586     stack_push((void *) im_data, IMRECT, im_free);
587     tv_init(imcalc_tv_get());
588     tv_init(imcal2_tv_get());
589     image_choice_reset();
590     fvector_free(data,imptrlz);
591     mixmodel_free(model_2d);
592 
593 }
594 
595 
596 void em_prob_class_plot_proc()
597 {
598    Imrect *im;
599 
600    if(prob_flag == 0)
601    {
602       error("Do <E-step> prior to probability calculation!", warning);
603       return;
604    }
605    im = prob_class(model,pim_array,em_str_name);
606    stack_push(im, IMRECT, im_free);
607    if (imcalc_tv_get()!=NULL) tv_init(imcalc_tv_get());
608    if (imcal2_tv_get()!=NULL) tv_init(imcal2_tv_get());
609    image_choice_reset();
610 }
611 
612 /*
613  *   Commented out - GAB 19 Feb 2003
614  *   Not completely deleted because it may be useful for testing.
615  *
616 static void prob_all_plot_proc()
617 {
618    if (model == NULL)
619    {
620       error("Input file is missing!", warning);
621       return;
622    }
623    if(prob_flag == 0)
624    {
625       error("Do <E-step> prior to probability calculation!", warning);
626       return;
627    }
628    prob_all(model,model->par_dens, pim_array);
629    imcalc_draw(imcalc_tv_get());
630    imcalc_draw(imcal2_tv_get());
631    image_choice_reset();
632 }
633  *
634  *  End of "commented out" for function prob_all_plot_proc - GAB 19 Feb 2003
635  */
636 
637 static void em_hist_proc()
638 {
639    em_hist_plot(em_get_model());
640 }
641 
642 void em_mean_proc()
643 {
644    mean_update(pim_array, tv_get_im_roi(seq_tv_get()), em_get_model());
645    prob_flag = 0;
646 }
647 
648 void em_cov_proc()
649 {
650    cov_update(pim_array, tv_get_im_roi(seq_tv_get()), em_get_model());
651    prob_flag = 0;
652 }
653 
654 void em_prior_proc()
655 {
656    priors_update(pim_array, tv_get_im_roi(seq_tv_get()), em_get_model());
657    prob_flag = 0;
658 }
659 
660 void em_model_image(void ***pim_arr)
661 {
662    Sequence *seq = seq_get_current();
663    Imrect *im_tot, *im_tot1,*im_tot2,*im_tot3,*im_tot_mean,*im;
664    int   lx, ux, ly, uy;
665    int loop,loop1,i,j,k;
666    int im_no;
667    Mixmodel *m = em_get_model();
668    Imrect *im_seq =  seq_image_get();
669    Imregion *roi = NULL;
670    void ***pim_mean = NULL;
671 
672    if (im_seq == NULL)
673    {
674       error("Image is missing!", warning);
675       return;
676    }
677    roi = im_seq->region;
678    if (roi == NULL)
679         return;
680    if (m == NULL)
681    {
682        error("Input file is missing!", warning);
683        return;
684    }
685    if (pim_arr == NULL)
686    {
687      error("Probability densities not calculated!", warning);
688        return;
689    }
690    lx = roi->lx;
691    ux = roi->ux;
692    ly = roi->ly;
693    uy = roi->uy;
694    im_no = get_seq_current_frame(seq);
695 
696    pim_mean = (void***)pvector_alloc(0,m->nmix);
697    for (i=0; i < model->nmix ; i++)
698       pim_mean[i] = (void *) im_alloc(uy-ly,ux-lx,roi,float_v);
699 
700    for(k = 0; k < m->nmix;k++)
701       for (i = ly; i < uy; ++i)
702          for (j = lx; j < ux; ++j)
703             im_put_pixf(m->vectors[k][im_no],(Imrect *)pim_mean[k],i,j);
704 
705    if(prob_flag == 0)
706    {
707       error("Do <E-step> prior to probability calculation!", warning);
708       return;
709    }
710    if( (im_tot=em_get_im_tot()) == NULL)
711    {
712         error("Normalisation image missing \n",warning);
713 /*
714       im_tot = prob_im_tot(m,pim_array,roi);
715 */
716        return;
717    }
718    im_tot3 = im_alloc(uy-ly, ux-lx, roi, float_v);
719    for(loop = 0; loop < m->nmix; loop++)
720    {
721       im_tot1 = im_alloc(uy-ly, ux-lx, roi, float_v);
722       for(loop1 = 0; loop1 < m->nmix; loop1++)
723       {
724          if (pim_array[loop][loop1] != NULL)
725          {
726             im = (Imrect *)pim_array[loop][loop1];
727             im_tot2 = im_sum(im_tot1,im);
728             im_free(im_tot1);
729             im_tot1=im_tot2;
730          }
731       }
732       im_tot2 = prob_im_norm(im_tot1,im_tot);
733       im_tot1 = imf_prod(im_tot2,(Imrect *)pim_mean[loop]);
734       im_tot_mean = im_sum(im_tot1,im_tot3);
735       im_free(im_tot3);
736       im_tot3 = im_tot_mean;
737       im_free(im_tot1);
738    }
739    stack_push(im_tot3, IMRECT, im_free);
740    imcalc_draw(imcalc_tv_get());
741    im_free(im_tot2);
742 
743    for (i=0; i < m->nmix ; i++)
744       im_free((Imrect *)pim_mean[i]);
745    pvector_free(pim_mean,m->nmix);
746 
747 }
748 
749 static void em_model_proc()
750 {
751    em_model_image(pim_array);
752 }
753 
754 /*************************************************************************/
755 /*           Estimation (E) step of EM algorithm                         */
756 /*    Probability densities are calculated for pure and mixture models   */
757 /*    For each model the density image is generated                      */
758 /*    The sum of all density at each pixel location is calculated        */
759 /*************************************************************************/
760 void em_estep_proc()
761 {
762 /*
763    Imregion *roi;
764    Imregion *roi = tv_get_im_roi(seq_tv_get());
765 */
766 
767    /* Flag set - density images to be calculated before probability images */
768    prob_flag = 1;
769    pim_array = em_estep(pim_array, model,0);
770 
771 }
772 
773 void em_estep_grad_proc()
774 {
775 /*
776    Imregion *roi;
777    Imregion *roi = tv_get_im_roi(seq_tv_get());
778 */
779 
780    /* Flag set - density images to be calculated before probability images */
781    prob_flag = 1;
782    pim_array = em_estep(pim_array, model,1);
783 
784 }
785 
786 
787 
788 /*********************************************************/
789 /* Loads data from parameter file into structure "model" */
790 /*********************************************************/
791 
792 static void input_proc(void)
793 {
794    extern Mixmodel *input_mixmodel(char *file_name);
795    int i,j;
796 
797    if (pim_array != NULL && model != 0)
798    {
799       for (i=0; i < model->nmix ; i++)
800          for (j=0; j< model->nmix; j++)
801             im_free(pim_array[i][j]);
802       parray_free(pim_array,0,0,model->nmix,model->nmix);
803    }
804 
805    model = input_mixmodel(file_name);
806 }
807 
808 /*********************************************************/
809 /* Stores parameters from structure "model" into a file  */
810 /*********************************************************/
811 
812 static void output_proc(void)
813 {
814     output_mixmodel(model,file_name);
815 }
816 
817 void nmr_fit_proc()
818 {
819     shistogram *peaks;
820     shistogram *nmr_scaled_fit(Imregion *roi, Mixmodel *m);
821     void graph_hfit(Tv * tv, shistogram * hist);
822     Tv *tv_graph = (Tv*)imcalc_graph_tv_get();
823     Imregion *roi = tv_get_im_roi(seq_tv_get());
824 
825     peaks = nmr_scaled_fit(roi, em_get_model());
826     graph_hfit(tv_graph, peaks);
827 }
828 
829 /* Runs the segmentation routine to convergence as tested in TINA Memo no. 2005-014 */
830 
831 void run_edgeseg_proc(void)
832 {
833    int i;
834 
835    for(i=0; i<2; i++)
836    {
837       em_estep_proc();
838       em_mean_proc();
839       em_prior_proc();
840    }
841    for(i=0; i<5; i++)
842    {
843       em_estep_proc();
844       em_mean_proc();
845       em_cov_proc();
846       em_prior_proc();
847    }
848    for(i=0; i<33; i++)
849    {
850       em_estep_grad_proc();
851       em_mean_proc();
852       em_cov_proc();
853       em_prior_proc();
854       kgrad_calc_proc();
855    }
856    em_estep_grad_proc();
857 }
858 
859 
860 void run_noedgeseg_proc(void)
861 {
862    int i;
863 
864    for(i=0; i<2; i++)
865    {
866       em_estep_proc();
867       em_mean_proc();
868       em_prior_proc();
869    }
870    for(i=0; i<38; i++)
871    {
872       em_estep_proc();
873       em_mean_proc();
874       em_cov_proc();
875       em_prior_proc();
876    }
877    em_estep_proc();
878 }
879 
880 
881 void  nmr_segment_tool(int x, int y)
882 {
883   static void *tool = NULL;
884 
885   if (tool)
886   {
887      tw_show_tool(tool);
888      return;
889   }
890 
891   tool = tw_tool("NMR Segment Tool", x, y);
892   pfile = (void *) tw_sglobal("Model File : ", file_name, 32);
893   tw_label("            ");
894   tw_help_button("nmr_segment_tool");
895   tw_newrow();
896   tw_button("input", input_proc, NULL);
897   tw_button("output", output_proc, NULL);
898   tw_button("ptest1", ptest_proc1, NULL);
899   tw_newrow();
900   tw_label("Multi-D fit: ");
901   tw_button("hfit-scale", nmr_fit_proc, NULL);
902   tw_button("E-step", em_estep_proc, NULL);
903   tw_button("E-step grad", em_estep_grad_proc, NULL);
904   tw_newrow();
905   tw_button("M-mean", em_mean_proc, NULL);
906   tw_button("M-cov", em_cov_proc, NULL);
907   tw_button("M-prior", em_prior_proc, NULL);
908   tw_button("M-k grad calc", kgrad_calc_proc, NULL);
909 
910   tw_newrow();
911   p_em_str_name = (void *) tw_sglobal("Tissue Type : ", em_str_name, 32);
912 
913   tw_newrow();
914   tw_label("EM plot :");
915   tw_button("1D hist", em_hist_proc, NULL);
916   tw_button("2D hist", ptest_proc, NULL);
917   tw_button("model", em_model_proc,NULL);
918   tw_button("prob", em_prob_class_plot_proc, NULL);
919 /*
920  *   Commented out - GAB 19 Feb 2003
921  *   Not completely deleted because it may be useful for testing.
922  *
923   tw_button("prob all", prob_all_plot_proc, NULL);
924 */
925 
926   tw_newrow();
927   tw_label("Muli-D Grad Calc:");
928   tw_button("slope image", mD_slope_proc, NULL);
929   tw_button("scat plot",   mD_grad_scat_proc, NULL);
930   tw_button("grad model",  md_grad_model_proc, NULL);
931   tw_newrow();
932   tw_button("Seg without grad", run_noedgeseg_proc, NULL);
933   tw_button("Seg with grad", run_edgeseg_proc, NULL);
934 
935   tw_end_tool();
936 }
937 

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