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

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

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