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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedSeg_nmrold_backend.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  *
  6  * Redistribution and use in source and binary forms, with or without modification,
  7  * are permitted provided that the following conditions are met:
  8  *
  9  *   . Redistributions of source code must retain the above copyright notice,
 10  *     this list of conditions and the following disclaimer.
 11  *
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation
 14  *     and/or other materials provided with the distribution.
 15  *
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this
 18  *     software without specific prior written permission.
 19  *
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :
 37  * Date    :
 38  * Version :
 39  * CVS Id  :
 40  *
 41  * Notes :
 42  *
 43 */
 44 
 45 #include "tlmedSeg_nmrold_backend.h"
 46 
 47 #if HAVE_CONFIG_H
 48   #include <config.h>
 49 #endif
 50 
 51 #include <stdio.h>
 52 #include <math.h>
 53 #include <tina/sys/sysDef.h>
 54 #include <tina/sys/sysPro.h>
 55 #include <tina/math/mathDef.h>
 56 #include <tina/math/mathPro.h>
 57 #include <tina/image/imgDef.h>
 58 #include <tina/image/imgPro.h>
 59 
 60 #include <tinatool/tlbase/tlbasePro.h>
 61 
 62 Prof1 *csfgmtri;
 63 Prof1 *csfgmgaus;
 64 Prof1 *csfgmlookup;
 65 Prof1 *gmwmtri;
 66 Prof1 *gmwmgaus;
 67 Prof1 *gmwmlookup;
 68 
 69 #define INTSTEP 0.2
 70 
 71 
 72 double dgauss(double a, double b, double c)
 73 {
 74    double temp;
 75    temp = (a-c)*(a-c)/(2.0*b*b);
 76    return(exp(-temp));
 77 }
 78 
 79 
 80 double lin_var(double var1,double var2,double x1,double x2, double x)
 81 {
 82     double fac,var;
 83     fac = (x-x1)/(x2-x1);
 84     if (fac<0.0) fac = 0.0;
 85     if (fac>1.0) fac = 1.0;
 86     var =  fabs(var2)*fac + fabs(var1)*(1.0-fac);
 87     return(var);
 88 }
 89 
 90 
 91 double nmr_gaus_tri(double sigma, double xs, double xl, double x)
 92 {
 93     double result=0.0;
 94     double step;
 95     double factor;
 96     double lastx;
 97 
 98     for (step = x-4.0*sigma; step< x+4.0*sigma; step+=INTSTEP*sigma)
 99     {
100         if ( ((step > xs) && (step <xl)) || ((step < xs) && (step > xl)))
101            result += (step-xs)*dgauss(x,sigma,step);
102         else if (fabs(step-xl)< INTSTEP*sigma)
103         {
104            factor = 1.0-fabs(step-xl)/(INTSTEP*sigma);
105            if (step > xl)
106            {
107              lastx = ((step-0.5*INTSTEP*sigma)+xl)/2.0;
108            }
109            else
110            {
111              lastx = ((step+0.5*INTSTEP*sigma)+xl)/2.0;
112            }
113            result += factor*(lastx-xs)*dgauss(x,sigma,lastx);
114         }
115     }
116     result *=INTSTEP/(sqrt(2.0*PI)*(xl-xs));
117     return(fabs(result));
118 }
119 
120 
121 float nmr_gm_f_prob(double *b, float pixel_val)
122 {
123 /* Works out the GM probabilities */
124 
125    float lowerlimit, upperlimit;
126    float new_pixel_val;
127    double var;
128    double gm_prob_gaus;
129    double gm_prob_csfgm;
130    double gm_prob_gmwm;
131 
132    /* Check if pixel_val is within the possible GM range i.e. */
133    /* 4 standard deviations either side of b[1] and b[8] */
134    lowerlimit = b[1]-4.0*b[2];
135    upperlimit = b[8]+4.0*b[9];
136 
137    if (pixel_val<(float)lowerlimit || pixel_val>(float)upperlimit)
138       new_pixel_val=0.0;
139    else
140    {
141       /* GM Probability due to Main Gaussian peak */
142       if (pixel_val>(b[4]+(4.0*b[5])) || pixel_val<(b[4]-(4.0*b[5])))
143          gm_prob_gaus=0.0;
144       else
145          gm_prob_gaus=b[3]*b[3]*dgauss(b[4],b[5],(double)pixel_val);
146 
147       /* GM Probability due to CSF GM Triangle */
148 
149       /* GM Probability due to GM WM Triangle */
150 
151       var = lin_var(b[5],b[2],b[4],b[1],(double)pixel_val);
152       gm_prob_csfgm = b[6]*b[6]*nmr_gaus_tri(var,b[1],b[4],(double)pixel_val);
153 
154       /* GM Probability due to GM WM Triangle */
155       var = lin_var(b[9],b[5],b[8],b[4],(double)pixel_val);
156       gm_prob_gmwm = b[10]*b[10]*nmr_gaus_tri(var,b[8],b[4],(double)pixel_val);
157 
158       new_pixel_val=gm_prob_gaus+gm_prob_csfgm+gm_prob_gmwm;
159    }
160 
161    return (new_pixel_val);
162 }
163 
164 
165 float nmr_wm_f_prob(double *b, float pixel_val)
166 {
167 /* Works out the WM probabilities */
168 
169    float lowerlimit;
170    float new_pixel_val;
171    double var;
172    double wm_prob_gaus;
173    double wm_prob_gmwm;
174 
175    /* Check if pixel_val is within the possible WM range i.e. */
176    /* greater than 4 standard deviations below of b[4] */
177    lowerlimit = b[4]-4.0*b[9];
178 
179    if (pixel_val<lowerlimit)
180       new_pixel_val=0.0;
181    else
182    {
183       /* WM Probability due to Main Gaussian peak */
184       if (pixel_val>(b[8]+(4.0*b[9])) || pixel_val<(b[8]-(4.0*b[9])))
185          wm_prob_gaus=0.0;
186       else
187          wm_prob_gaus=b[7]*b[7]*dgauss(b[8],b[9],(double)pixel_val);
188 
189       /* WM Probability due to GM WM Triangle */
190       var = lin_var(b[5],b[9],b[4],b[8], (double)pixel_val);
191       wm_prob_gmwm=b[10]*b[10]*nmr_gaus_tri(var,b[4],b[8], (double)pixel_val);
192 
193       new_pixel_val=wm_prob_gaus+wm_prob_gmwm;
194    }
195 
196    return (new_pixel_val);
197 }
198 
199 
200 float nmr_csf_f_prob(double *b, float pixel_val)
201 {
202 /* Works out the CSF probabilities */
203 
204    float upperlimit;
205    float new_pixel_val;
206    double var;
207    double csf_prob_gaus;
208    double csf_prob_csfgm;
209 
210    /* Check if pixel_val is within the possible CSF range i.e. */
211    /* less than 4 standard deviations above b[4] */
212    upperlimit = b[4]+4.0*b[2];
213 
214    if (pixel_val>upperlimit)
215       new_pixel_val=0.0;
216    else
217    {
218       /* CSF Probability due to Main Gaussian peak */
219       if (pixel_val>(b[1]+(4.0*b[2])) || pixel_val<(b[1]-(4.0*b[2])))
220          csf_prob_gaus=0.0;
221       else
222          csf_prob_gaus=b[0]*b[0]*dgauss(b[1],b[2],(double)pixel_val);
223 
224       /* CSF Probability due to CSF GM Triangle */
225       var = lin_var(b[5],b[2],b[4],b[1],(double)pixel_val);
226       csf_prob_csfgm=b[6]*b[6]*(nmr_gaus_tri(var,b[4],b[1],(double)pixel_val));
227 
228       new_pixel_val=csf_prob_gaus+csf_prob_csfgm;
229    }
230 
231    return (new_pixel_val);
232 }
233 
234 
235 static Imrect *nmr_hfit_gm_f(double *b, Imrect * im1)
236 {
237 /* Processing for GM segmentation */
238 
239     Imrect *im2;
240     Imregion *roi;
241     float  *row1, *row2;
242     int     lx, ux, ly, uy;
243     int     i, j;
244  
245     if (im1 == NULL)
246         return (NULL);
247 
248     roi = im1->region;
249     if (roi == NULL)
250         return (NULL);
251 
252     lx = roi->lx;
253     ux = roi->ux;
254     ly = roi->ly;
255     uy = roi->uy;
256  
257     im2 = im_alloc(im1->height, im1->width, roi, float_v);
258 
259 
260     row1 = fvector_alloc(lx, ux);
261     row2 = fvector_alloc(lx, ux);
262 
263     for (i = ly; i < uy; ++i)
264     {
265         im_get_rowf(row1, im1, i, lx, ux);
266 
267         for (j = lx; j < ux; ++j)
268                 row2[j] = nmr_gm_f_prob(b, row1[j]);
269 
270         im_put_rowf(row2, im2, i, lx, ux);
271     }
272 
273     fvector_free(row1, lx);
274     fvector_free(row2, lx);
275 
276     return (im2);
277 }
278 
279 
280 static Imrect *nmr_hfit_wm_f(double *b, Imrect * im1)
281 /* Processing for WM segmentation */
282 {
283     Imrect *im3;
284     Imregion *roi;
285     float  *row1, *row2;
286     int     lx, ux, ly, uy;
287     int     i, j;
288  
289     if (im1 == NULL)
290         return (NULL);
291 
292 
293     roi = im1->region;
294     if (roi == NULL)
295         return (NULL);
296 
297     lx = roi->lx;
298     ux = roi->ux;
299     ly = roi->ly;
300     uy = roi->uy;
301  
302 
303     im3 = im_alloc(im1->height, im1->width, roi, float_v);
304 
305     row1 = fvector_alloc(lx, ux);
306     row2 = fvector_alloc(lx, ux);
307 
308     for (i = ly; i < uy; ++i)
309     {
310         im_get_rowf(row1, im1, i, lx, ux);
311 
312         for (j = lx; j < ux; ++j)
313                 row2[j] = nmr_wm_f_prob(b, row1[j]);
314 
315         im_put_rowf(row2, im3, i, lx, ux);
316     }
317     fvector_free(row1, lx);
318     fvector_free(row2, lx);
319 
320     return (im3);
321 }
322 
323 
324 static Imrect *nmr_hfit_csf_f(double *b, Imrect * im1)
325 {
326 /* Processing for CSF segmentation */
327 
328     Imrect *im4;
329     Imregion *roi;
330     float  *row1, *row2;
331     int     lx, ux, ly, uy;
332     int     i, j;
333  
334 
335     if (im1 == NULL)
336         return (NULL);
337 
338     roi = im1->region;
339     if (roi == NULL)
340         return (NULL);
341 
342     lx = roi->lx;
343     ux = roi->ux;
344     ly = roi->ly;
345     uy = roi->uy;
346  
347 
348     im4 = im_alloc(im1->height, im1->width, roi, float_v);
349 
350 
351     row1 = fvector_alloc(lx, ux);
352     row2 = fvector_alloc(lx, ux);
353 
354     for (i = ly; i < uy; ++i)
355     {
356         /*get the row from im1*/
357         im_get_rowf(row1, im1, i, lx, ux);
358 
359         for (j = lx; j < ux; ++j)
360                 row2[j] = nmr_csf_f_prob(b, row1[j]);
361 
362         im_put_rowf(row2, im4, i, lx, ux);
363     }
364 
365     /*free the rows*/ 
366     fvector_free(row1, lx);
367     fvector_free(row2, lx);
368 
369     return (im4);
370 }
371 
372 
373 static Imrect *nmr_hfit_prob_total(Imrect *im2, Imrect *im3, Imrect *im4)
374 {
375 /* Totalling all the image probabilities - for normalisation */
376 
377    Imrect *im5,*im6;
378 
379    im5 = im_sum(im2,im3);
380    im6 = im_sum(im4,im5);
381    im_free(im5);
382 
383    return (im6);
384 }
385 
386 
387 static Imrect *nmr_hfit_normed(Imrect *im_orig, Imrect *im_tot)
388 {
389 /* Normalisation of each segmented image */
390 
391    Imrect *im_norm;
392    float threshold=1.0;
393    float rep_val=1.0;
394 
395    im_norm = im_div(im_orig,im_tot,threshold,rep_val);
396 
397    return (im_norm);
398 }
399 
400 
401 void gaus_and_blur_tri_func(double *b, int select)
402 {
403     Imrect *im;
404     Imrect *im2;
405     Imrect *im3;
406     Imrect *im4;
407     Imrect *im5;
408     int type;
409     Tv *imcalc_graph_tv_get();
410 
411  
412     if (stack_check_types(IMRECT, NULL) == false)
413     {
414         error("gaus_and_blur_tri_func: wrong type on stack", warning);
415         return;
416     }
417 
418     im = (Imrect *) stack_inspect(&type);
419 
420 
421     switch (im->vtype)
422     {
423        case complex_v:                  /*For complex images */
424           format("gaus_and_blur_tri_func: wrong type on stack", warning);
425           break;
426        default:                         /*For all other images */
427           im2 = nmr_hfit_gm_f(b, im);   /*Work out the grey matter */
428           im3 = nmr_hfit_wm_f(b, im);   /*Work out the white matter */
429           im4 = nmr_hfit_csf_f(b, im);  /*Work out the CSF */
430 
431           im5 = nmr_hfit_prob_total(im2,im3,im4); /*Work out the total probs */
432           switch (select)
433           {
434              case 0:
435                im = nmr_hfit_normed(im4,im5); /*Normalised CSF image */
436              break;
437              case 1:
438                im = nmr_hfit_normed(im2,im5); /*Normalised GM image */
439              break;
440              case 2:
441                im = nmr_hfit_normed(im3,im5); /*Normalised WM image */
442              break;
443           }
444           im_free(im2);
445           im_free(im3);
446           im_free(im4);
447           im_free(im5);
448 /*
449           stack_push((void *) im4, IMRECT, im_free);
450           stack_push((void *) im3, IMRECT, im_free);
451           stack_push((void *) im2, IMRECT, im_free);
452 */
453           stack_push((void *) im, IMRECT, im_free);
454     }
455 
456 
457     imcalc_draw(imcalc_tv_get());
458     imcalc_draw(imcal2_tv_get());
459 
460     image_choice_reset();
461 }
462 

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