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

Linux Cross Reference
Tina4/src/tools/nmr-segment/gaus_blur_tri_seg.c

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

  1 #include <stdio.h>
  2 #include <sys/param.h>
  3 #include <string.h>
  4 #include <tina/sys.h>
  5 #include <tina/sysfuncs.h>
  6 #include <tina/math.h>
  7 #include <tina/mathfuncs.h>
  8 #include <tina/vision.h>
  9 #include <tina/visionfuncs.h>
 10 #include <tina/file.h>
 11 #include <tina/filefuncs.h>
 12 #include <tina/file_gen.h>
 13 #include <tina/tv.h>
 14 #include <tina/tvfuncs.h>
 15 #include <tina/draw.h>
 16 #include <tina/drawfuncs.h>
 17 #include <tina/toolsfuncs.h>
 18 #include <tina/tw_Xfuncs.h>
 19 #include <tina/hist_funcs.h>
 20 
 21 Prof1 *csfgmtri;
 22 Prof1 *csfgmgaus;
 23 Prof1 *csfgmlookup;
 24 Prof1 *gmwmtri;
 25 Prof1 *gmwmgaus;
 26 Prof1 *gmwmlookup;
 27 extern double dgauss(double a, double b, double c);
 28 extern double norm_integ(double sigma,double mu,double y);
 29 extern void graph_hfit(Tv *tv, shistogram * hist);
 30 extern void csfgmtrigen (double start, double end, double height, int test_flag);
 31 extern void gmwmtrigen (double start, double end, double height, int test_flag);
 32 extern Prof1 *gausgen (double stddev, int test_flag);
 33 extern Prof1 *convolve (Prof1 *arr1, Prof1 *arr2, int test_flag);
 34 extern void graph_hfit(Tv *tv, shistogram * hist);
 35 
 36 #define INTSTEP 0.2
 37 
 38 double lin_var(double var1,double var2,double x1,double x2, double x)
 39 {
 40     double fac,var;
 41     fac = (x-x1)/(x2-x1);
 42     if (fac<0.0) fac = 0.0;
 43     if (fac>1.0) fac = 1.0;
 44     var =  fabs(var2)*fac + fabs(var1)*(1.0-fac);
 45     return(var);
 46 }
 47 
 48 double nmr_gaus_tri(double sigma, double xs, double xl, double x)
 49 {
 50     double result=0.0;
 51     double step;
 52     double factor;
 53     double lastx;
 54 
 55     for (step = x-4.0*sigma; step< x+4.0*sigma; step+=INTSTEP*sigma)
 56     {
 57         if ( ((step > xs) && (step <xl)) || ((step < xs) && (step > xl)))
 58            result += (step-xs)*dgauss(x,sigma,step);
 59         else if (fabs(step-xl)< INTSTEP*sigma)
 60         {
 61            factor = 1.0-fabs(step-xl)/(INTSTEP*sigma);
 62            if (step > xl)
 63            {
 64              lastx = ((step-0.5*INTSTEP*sigma)+xl)/2.0;
 65            }
 66            else
 67            {
 68              lastx = ((step+0.5*INTSTEP*sigma)+xl)/2.0;
 69            }
 70            result += factor*(lastx-xs)*dgauss(x,sigma,lastx);
 71         }
 72     }
 73     result *=INTSTEP/(sqrt(2.0*PI)*(xl-xs));
 74 
 75 
 76 /*
 77     result = sigma*(dgauss(x,sigma,xs) - dgauss(x,sigma,xl))
 78             /(sqrt(2.0*PI)*fabs(xl-xs));
 79     result += fabs((x - xs)*(norm_integ(sigma,xl,x) 
 80                       - norm_integ(sigma,xs,x))/(xl-xs));
 81 */
 82     return(fabs(result));
 83 }
 84 
 85 float nmr_gm_f_prob(double *b, float pixel_val)
 86 {
 87 /* Works out the GM probabilities */
 88 
 89    float lowerlimit, upperlimit;
 90    float new_pixel_val;
 91    double var;
 92    double gm_prob_gaus;
 93    double gm_prob_csfgm;
 94    double gm_prob_gmwm;
 95 
 96    /* Check if pixel_val is within the possible GM range i.e. */
 97    /* 4 standard deviations either side of b[1] and b[8] */
 98    lowerlimit = b[1]-4.0*b[2];
 99    upperlimit = b[8]+4.0*b[9];
100 
101    if (pixel_val<(float)lowerlimit || pixel_val>(float)upperlimit)
102       new_pixel_val=0.0;
103    else
104    {
105       /* GM Probability due to Main Gaussian peak */
106       if (pixel_val>(b[4]+(4.0*b[5])) || pixel_val<(b[4]-(4.0*b[5])))
107          gm_prob_gaus=0.0;
108       else
109          gm_prob_gaus=b[3]*b[3]*dgauss(b[4],b[5],(double)pixel_val);
110 
111       /* GM Probability due to CSF GM Triangle */
112 
113       /* GM Probability due to GM WM Triangle */
114 
115       var = lin_var(b[5],b[2],b[4],b[1],(double)pixel_val);
116       gm_prob_csfgm = b[6]*b[6]*nmr_gaus_tri(var,b[1],b[4],(double)pixel_val);
117 
118       /* GM Probability due to GM WM Triangle */
119       var = lin_var(b[9],b[5],b[8],b[4],(double)pixel_val);
120       gm_prob_gmwm = b[10]*b[10]*nmr_gaus_tri(var,b[8],b[4],(double)pixel_val);
121 
122       new_pixel_val=gm_prob_gaus+gm_prob_csfgm+gm_prob_gmwm;
123    }
124 
125    return (new_pixel_val);
126 }
127 
128 
129 float nmr_wm_f_prob(double *b, float pixel_val)
130 {
131 /* Works out the WM probabilities */
132 
133    float lowerlimit;
134    float new_pixel_val;
135    double var;
136    double wm_prob_gaus;
137    double wm_prob_gmwm;
138 
139    /* Check if pixel_val is within the possible WM range i.e. */
140    /* greater than 4 standard deviations below of b[4] */
141    lowerlimit = b[4]-4.0*b[9];
142 
143    if (pixel_val<lowerlimit)
144       new_pixel_val=0.0;
145    else
146    {
147       /* WM Probability due to Main Gaussian peak */
148       if (pixel_val>(b[8]+(4.0*b[9])) || pixel_val<(b[8]-(4.0*b[9])))
149          wm_prob_gaus=0.0;
150       else
151          wm_prob_gaus=b[7]*b[7]*dgauss(b[8],b[9],(double)pixel_val);
152 
153       /* WM Probability due to GM WM Triangle */
154       var = lin_var(b[5],b[9],b[4],b[8], (double)pixel_val);
155       wm_prob_gmwm=b[10]*b[10]*nmr_gaus_tri(var,b[4],b[8], (double)pixel_val);
156 
157       new_pixel_val=wm_prob_gaus+wm_prob_gmwm;
158    }
159 
160    return (new_pixel_val);
161 }
162 
163 float nmr_csf_f_prob(double *b, float pixel_val)
164 {
165 /* Works out the CSF probabilities */
166 
167    float upperlimit;
168    float new_pixel_val;
169    double var;
170    double csf_prob_gaus;
171    double csf_prob_csfgm;
172 
173    /* Check if pixel_val is within the possible CSF range i.e. */
174    /* less than 4 standard deviations above b[4] */
175    upperlimit = b[4]+4.0*b[2];
176 
177    if (pixel_val>upperlimit)
178       new_pixel_val=0.0;
179    else
180    {
181       /* CSF Probability due to Main Gaussian peak */
182       if (pixel_val>(b[1]+(4.0*b[2])) || pixel_val<(b[1]-(4.0*b[2])))
183          csf_prob_gaus=0.0;
184       else
185          csf_prob_gaus=b[0]*b[0]*dgauss(b[1],b[2],(double)pixel_val);
186 
187       /* CSF Probability due to CSF GM Triangle */
188       var = lin_var(b[5],b[2],b[4],b[1],(double)pixel_val);
189       csf_prob_csfgm=b[6]*b[6]*(nmr_gaus_tri(var,b[4],b[1],(double)pixel_val));
190 
191       new_pixel_val=csf_prob_gaus+csf_prob_csfgm;
192    }
193 
194    return (new_pixel_val);
195 }
196 
197 static Imrect *nmr_hfit_gm_f(double *b, Imrect * im1)
198 {
199 /* Processing for GM segmentation */
200 
201     Imrect *im2;
202     Imregion *roi;
203     float  *row1, *row2;
204     int     lx, ux, ly, uy;
205     int     i, j;
206  
207     if (im1 == NULL)
208         return (NULL);
209 
210     roi = im1->region;
211     if (roi == NULL)
212         return (NULL);
213 
214     lx = roi->lx;
215     ux = roi->ux;
216     ly = roi->ly;
217     uy = roi->uy;
218  
219     im2 = im_alloc(im1->height, im1->width, roi, float_v);
220 
221 
222     row1 = fvector_alloc(lx, ux);
223     row2 = fvector_alloc(lx, ux);
224 
225     for (i = ly; i < uy; ++i)
226     {
227         im_get_rowf(row1, im1, i, lx, ux);
228 
229         for (j = lx; j < ux; ++j)
230                 row2[j] = nmr_gm_f_prob(b, row1[j]);
231 
232         im_put_rowf(row2, im2, i, lx, ux);
233     }
234 
235     fvector_free((void *) row1, lx);
236     fvector_free((void *) row2, lx);
237 
238     return (im2);
239 }
240 
241 
242 static Imrect *nmr_hfit_wm_f(double *b, Imrect * im1)
243 /* Processing for WM segmentation */
244 {
245     Imrect *im3;
246     Imregion *roi;
247     float  *row1, *row2;
248     int     lx, ux, ly, uy;
249     int     i, j;
250  
251     if (im1 == NULL)
252         return (NULL);
253 
254 
255     roi = im1->region;
256     if (roi == NULL)
257         return (NULL);
258 
259     lx = roi->lx;
260     ux = roi->ux;
261     ly = roi->ly;
262     uy = roi->uy;
263  
264 
265     im3 = im_alloc(im1->height, im1->width, roi, float_v);
266 
267     row1 = fvector_alloc(lx, ux);
268     row2 = fvector_alloc(lx, ux);
269 
270     for (i = ly; i < uy; ++i)
271     {
272         im_get_rowf(row1, im1, i, lx, ux);
273 
274         for (j = lx; j < ux; ++j)
275                 row2[j] = nmr_wm_f_prob(b, row1[j]);
276 
277         im_put_rowf(row2, im3, i, lx, ux);
278     }
279     fvector_free((void *) row1, lx);
280     fvector_free((void *) row2, lx);
281 
282     return (im3);
283 }
284 
285 static Imrect *nmr_hfit_csf_f(double *b, Imrect * im1)
286 {
287 /* Processing for CSF segmentation */
288 
289     Imrect *im4;
290     Imregion *roi;
291     float  *row1, *row2;
292     int     lx, ux, ly, uy;
293     int     i, j;
294  
295 
296     if (im1 == NULL)
297         return (NULL);
298 
299     roi = im1->region;
300     if (roi == NULL)
301         return (NULL);
302 
303     lx = roi->lx;
304     ux = roi->ux;
305     ly = roi->ly;
306     uy = roi->uy;
307  
308 
309     im4 = im_alloc(im1->height, im1->width, roi, float_v);
310 
311 
312     row1 = fvector_alloc(lx, ux);
313     row2 = fvector_alloc(lx, ux);
314 
315     for (i = ly; i < uy; ++i)
316     {
317         /*get the row from im1*/
318         im_get_rowf(row1, im1, i, lx, ux);
319 
320         for (j = lx; j < ux; ++j)
321                 row2[j] = nmr_csf_f_prob(b, row1[j]);
322 
323         im_put_rowf(row2, im4, i, lx, ux);
324     }
325 
326     /*free the rows*/ 
327     fvector_free((void *) row1, lx);
328     fvector_free((void *) row2, lx);
329 
330     return (im4);
331 }
332 
333 
334 static Imrect *nmr_hfit_prob_total(Imrect *im2, Imrect *im3, Imrect *im4)
335 {
336 /* Totalling all the image probabilities - for normalisation */
337 
338    Imrect *im5,*im6;
339 
340    im5 = im_sum(im2,im3);
341    im6 = im_sum(im4,im5);
342    im_free(im5);
343 
344    return (im6);
345 }
346 
347 
348 static Imrect *nmr_hfit_normed(Imrect *im_orig, Imrect *im_tot)
349 {
350 /* Normalisation of each segmented image */
351 
352    Imrect *im_norm;
353    float threshold=1.0;
354    float rep_val=1.0;
355 
356    im_norm = im_div(im_orig,im_tot,threshold,rep_val);
357 
358    return (im_norm);
359 }
360 
361 
362 void gaus_and_blur_tri_func(double *b, int select)
363 {
364     Imrect *im;
365     Imrect *im2;
366     Imrect *im3;
367     Imrect *im4;
368     Imrect *im5;
369     int type;
370     Tv *imcalc_graph_tv_get();
371 
372  
373     if (stack_check_types(IMRECT, NULL) == false)
374     {
375         error("gaus_and_blur_tri_func: wrong type on stack", warning);
376         return;
377     }
378 
379     im = (Imrect *) stack_inspect(&type);
380 
381 
382     switch (im->vtype)
383     {
384        case complex_v:                  /*For complex images */
385           format("gaus_and_blur_tri_func: wrong type on stack", warning);
386           break;
387        default:                         /*For all other images */
388           im2 = nmr_hfit_gm_f(b, im);   /*Work out the grey matter */
389           im3 = nmr_hfit_wm_f(b, im);   /*Work out the white matter */
390           im4 = nmr_hfit_csf_f(b, im);  /*Work out the CSF */
391 
392           im5 = nmr_hfit_prob_total(im2,im3,im4); /*Work out the total probs */
393           switch (select)
394           {
395              case 0:
396                im = nmr_hfit_normed(im4,im5); /*Normalised CSF image */
397              break;
398              case 1:
399                im = nmr_hfit_normed(im2,im5); /*Normalised GM image */
400              break;
401              case 2:
402                im = nmr_hfit_normed(im3,im5); /*Normalised WM image */
403              break;
404           }
405           im_free(im2);
406           im_free(im3);
407           im_free(im4);
408           im_free(im5);
409 /*
410           stack_push((void *) im4, IMRECT, im_free);
411           stack_push((void *) im3, IMRECT, im_free);
412           stack_push((void *) im2, IMRECT, im_free);
413 */
414           stack_push((void *) im, IMRECT, im_free);
415     }
416 
417 
418     imcalc_draw(imcalc_tv_get());
419     imcalc_draw(imcal2_tv_get());
420 
421     image_choice_reset();
422 }
423 

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