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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedSeg_nmrold.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.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 #include <tinatool/tlmedical/tlmedSeg_hist.h>
 60 
 61 
 62 /* interactive global variables for access via dialog box */
 63 double *a;
 64 /*
 65 void *pcsf_sigma;
 66 void *pgrey_sigma;
 67 void *pwhite_sigma;
 68 void *pcsf_peak;
 69 void *pgrey_peak;
 70 void *pwhite_peak;
 71 void *pcsf_peak1;
 72 void *pgrey_peak1;
 73 void *pwhite_peak1;
 74 void *pcsf_peak2;
 75 void *pgrey_peak2;
 76 void *pwhite_peak2;
 77 */
 78 double csf_sigma;
 79 double grey_sigma;
 80 double white_sigma;
 81 double csf_peak;
 82 double grey_peak;
 83 double white_peak;
 84 double csf_min;
 85 double white_max;
 86 double p1c;
 87 double p1g;
 88 double p1w;
 89 double p2c;
 90 double p2g;
 91 double p2w;
 92 int seg_im;
 93 
 94 shistogram *peaks;
 95 
 96 double glob_sigma; /* Global sigma for the smoothing kernel for GNC */
 97 double glob_acc;   /* Global acc for the smoothing kernal for GNC */
 98 int glob_smooth_flag; /* Whether to smooth or not in hist_fit_pab */
 99 
100 
101 double nmr_hfit_old(int n_par, double *a, float x)
102 {
103    double temp = 0.0;
104    float nmr_wm_f_prob(double *b, float pixel_val);
105    float nmr_gm_f_prob(double *b, float pixel_val);
106    float nmr_csf_f_prob(double *b, float pixel_val);
107 
108    a[2] = a[9];
109    if (a[1] > a[4])
110    {
111       a[3] = 0.0;
112       a[4] = a[1];
113    }
114    if (a[4] > a[8])
115    {
116       a[7] = 0.0;
117       a[8] = a[4];
118    }
119    temp +=  nmr_wm_f_prob(a, x)
120          + nmr_csf_f_prob(a, x)
121          + nmr_gm_f_prob(a, x);
122    return(temp);
123 }
124 
125 
126 static double pixchisq2_old(int n_par, double *a)
127 {
128   double chisq = 0.0, temp;
129   int     i;
130   float x;
131 
132   for (i=0;i<peaks->xbins;i++)
133     {
134       x = peaks->xmin+(i+0.5)*peaks->xincr;
135       temp = nmr_hfit_old(n_par, a , x);
136       chisq += (peaks->array[0][i]-temp)*(peaks->array[0][i]-temp)/(temp+1.0);
137     }
138   return(chisq);
139 }
140 
141 
142 void hist_fit_old(shistogram *hist)
143 {
144     int n = 11;
145     double c_test1 = 0.00001;
146     double grey_vol, white_vol, csf_vol;
147     double grey_white, csf_grey;
148     double *lambda;
149 
150     if (a!=NULL) rfree(a);
151     a = (double *)ralloc(n*sizeof(double));
152     lambda = (double *)ralloc(n*sizeof(double));
153     peaks = hist;
154 
155     a[1] = csf_peak;
156     a[4] = grey_peak;
157     a[8] = white_peak;
158     a[6] = hfill1(hist,(float)((a[1]+a[4])/2.0), 0.0);
159     a[6] = sqrt(a[6]);
160     a[0] = hfill1(hist, (float)a[1], 0.0)-a[6]*a[6];
161     if (a[0]>0.0) a[0] = sqrt(a[0]);
162     else a[0] = 0.0;
163     a[3] = hfill1(hist, (float)a[4], 0.0)-a[6]*a[6];
164     if(a[3]>0.0) a[3] = sqrt(a[3]);
165     else a[3] = 0.0;
166     a[7] = hfill1(hist, (float)a[8], 0.0);
167     a[7] = sqrt(a[7]);
168     a[2] = csf_sigma;
169     a[5] = grey_sigma;
170     a[9] = white_sigma;
171     a[10] = 2.0;
172     lambda[0] = 10.0;
173     lambda[1] = 4.0;
174     lambda[2] = 1.0;
175     lambda[3] = 10.0;
176     lambda[4] = 4.0;
177     lambda[5] = 1.0;
178     lambda[6] = 10.0;
179     lambda[7] = 10.0;
180     lambda[8] = 4.0;
181     lambda[9] = 1.0;
182     lambda[10] = 5.0;
183 
184     simplexmin2(n, a, lambda, pixchisq2_old, NULL, c_test1, (void (*) (char *)) format);
185     hsuper(hist, n, nmr_hfit_old, a);
186 
187     csf_peak = a[1];
188     grey_peak = a[4];
189     white_peak = a[8];
190     csf_sigma = a[2];
191     grey_sigma = a[5];
192     white_sigma = a[9];
193     csf_vol = sqrt(2.0*3.141592)*a[0]*a[0]*a[2]/hist->xincr;
194     csf_grey  = (a[4]-a[1])*a[6]*a[6]/hist->xincr;
195     grey_vol = sqrt(2.0*3.141592)*a[3]*a[3]*a[5]/hist->xincr;
196     grey_white = (a[8]-a[4])*a[10]*a[10]/hist->xincr;
197     white_vol = sqrt(2.0*3.141592)*a[7]*a[7]*a[9]/hist->xincr;
198     format(" csf   = %3.2f + %3.2f = %3.2f \n", csf_vol,csf_grey/2.0, csf_vol+csf_grey/2.0);
199     format(" grey  = %3.2f + %3.2f + %3.2f = %3.2f \n", grey_vol, csf_grey/2.0,
200            grey_white/2.0, grey_vol + csf_grey/2.0+ grey_white/2.0);
201     format(" white = %3.2f + %3.2f = %3.2f \n", white_vol, grey_white/2.0, white_vol+grey_white/2.0);
202     rfree(lambda);
203     return;
204 }
205 
206 
207 double nmr_hfit_pab(int n_par, double *a, float x)
208 {
209 /* Function to fit just the CSF peak and partial volume linear component for a T2 IR brain MRI */
210 /* Assume a[0] = amplitude of csf peak
211           a[1] = mean of csf peak
212           a[2] = s.d. of csf peak
213           a[3] = s.d. of partial volume peak
214           a[4] = amplitude of partial volume peak
215 */
216         double dgauss(double a, double b, double c);
217         double erf(double x);
218         double b, b1, c, c1, z, tot;
219 
220         z = (a[1]-x)/(sqrt(2)*a[2]);
221         b = dgauss(a[1], a[2], (double)x);
222         b1 = a[0]*a[0]*b;
223         c = erf(z);
224         c1 = a[4]*a[4]*(0.5-0.5*c);
225         tot = b1+c1;
226 
227         return tot;
228 }
229 
230 
231 float *make_model_vec_pab(int n_par, double *a)
232 {
233 
234 /* Having to make this all float since the profile smoothing functions work with floats: THIS IS A BAD THING */
235 
236   float *model_vec=NULL;
237   int     i;
238   float x;
239   Prof1 *prof=NULL;
240 
241     model_vec = fvector_alloc(0, peaks->xbins);
242 
243     for (i=0;i<peaks->xbins;i++)
244     {
245       x = peaks->xmin+(i+0.5)*peaks->xincr;
246       model_vec[i] = (float)nmr_hfit_pab(n_par, a , x);
247     }
248 
249     if(glob_sigma!=0.0)
250     {
251         prof = prof_gauss(glob_sigma, glob_acc);
252         smooth_1d_sym(model_vec, 0, peaks->xbins, prof);
253     }
254     prof1_free(prof);
255     return model_vec;
256 }
257 
258 
259 double pixchisq2_pab(int n_par, double *a)
260 {
261         double chisq = 0.0, temp;
262         int i;
263         float x, *model_vec=NULL;
264 
265         if(glob_smooth_flag==1)
266         {
267                 model_vec = make_model_vec_pab(n_par, a);
268                 for (i=0;i<peaks->xbins;i++)
269                 {
270                         chisq += (peaks->array[0][i]-model_vec[i])*(peaks->array[0][i]-model_vec[i])/(model_vec[i]+1.0);
271                 }
272                 fvector_free(model_vec, 0);
273         }
274         else
275         {
276                 for (i=0;i<peaks->xbins;i++)
277                 {
278                         x = peaks->xmin+(i+0.5)*peaks->xincr;
279                         temp = nmr_hfit_pab(n_par, a , x);
280                         chisq += (peaks->array[0][i]-temp)*(peaks->array[0][i]-temp)/(temp+1.0);
281                 }
282         }
283         return(chisq);
284 }
285 
286 
287 void hist_fit_pab(shistogram *hist, int smooth_flag, double smooth_sigma )
288 {
289     int n = 5, i;
290     double c_test1 = 0.00001;
291     double *lambda;
292     float *temp=NULL;
293     double acc;
294     Prof1 *prof=NULL;
295 
296     if (a!=NULL) rfree(a);
297     a = (double *)ralloc(n*sizeof(double));
298     lambda = (double *)ralloc(n*sizeof(double));
299     peaks = hist;
300 
301     a[1] = csf_peak;
302     a[0] = hfill1(hist, (float)a[1], 0.0);
303     a[0] = sqrt(a[0]);
304     a[2] = csf_sigma;
305     a[3] = csf_sigma;
306     a[4] = hfill1(hist, (float)(white_max-1), 0.0);  /* The minus 1 is to make sure you stay in the histogram range */
307     a[4] = sqrt(a[4]);
308 
309     lambda[0] = 10.0;
310     lambda[1] = 4.0;
311     lambda[2] = 1.0;
312     lambda[3] = 1.0;
313     lambda[4] = 10.0;
314 
315     glob_smooth_flag = smooth_flag;
316 
317     if(smooth_flag==1)
318     {
319         temp = fvector_alloc(0, peaks->xbins);
320         for(i=0; i<peaks->xbins; i++)
321         {
322             temp[i] = peaks->array[0][i];
323         }
324 
325         acc = 0.1;
326         glob_sigma = smooth_sigma;
327         glob_acc = acc;
328 
329         prof = prof_gauss(smooth_sigma, acc);
330         smooth_1d_sym(peaks->array[0], 0, peaks->xbins, prof);
331         prof1_free(prof);
332     }
333 
334     simplexmin2(n, a, lambda, pixchisq2_pab, NULL, c_test1, (void (*) (char *)) format);
335 
336     if(smooth_flag==1)
337     {
338 /*
339         for(i=0; i<peaks->xbins; i++)
340         {
341             peaks->array[0][i] = temp[i];
342         }
343 */
344         fvector_free(temp, 0);
345     }
346 
347     hsuper(hist, n, nmr_hfit_pab, a);
348 
349     csf_peak = a[1];
350     csf_sigma = a[2];
351     rfree(lambda);
352     return;
353 }
354 

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