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

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

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

  1 /*
  2   nmr segmentation (static image)
  3 
  4   n.thacker
  5 */
  6 
  7 #include <stdio.h>
  8 #include <sys/param.h>
  9 #include <string.h>
 10 #include <tina/sys.h>
 11 #include <tina/sysfuncs.h>
 12 #include <tina/math.h>
 13 #include <tina/mathfuncs.h>
 14 #include <tina/vision.h>
 15 #include <tina/visionfuncs.h>
 16 #include <tina/file.h>
 17 #include <tina/filefuncs.h>
 18 #include <tina/file_gen.h>
 19 #include <tina/tv.h>
 20 #include <tina/tvfuncs.h>
 21 #include <tina/draw.h>
 22 #include <tina/drawfuncs.h>
 23 #include <tina/toolsfuncs.h>
 24 #include <tina/tw_Xfuncs.h>
 25 #include <tina/hist_funcs.h>
 26 
 27 
 28 static int param_num = 1;
 29 
 30 /* interactive global variables for access via dialog box */
 31 extern  double *a;
 32 extern  void *pcsf_sigma;
 33 extern  void *pgrey_sigma;
 34 extern  void *pwhite_sigma;
 35 extern  void *pcsf_peak;
 36 extern  void *pgrey_peak;
 37 extern  void *pwhite_peak;
 38 extern  void *pcsf_peak1;
 39 extern  void *pgrey_peak1;
 40 extern  void *pwhite_peak1;
 41 extern  void *pcsf_peak2;
 42 extern  void *pgrey_peak2;
 43 extern  void *pwhite_peak2;
 44 extern double csf_sigma;
 45 extern double grey_sigma;
 46 extern double white_sigma;
 47 extern double csf_peak;
 48 extern double grey_peak;
 49 extern double white_peak;
 50 extern double csf_min;
 51 extern double white_max;
 52 extern double p1c;
 53 extern double p1g;
 54 extern double p1w;
 55 extern double p2c;
 56 extern double p2g;
 57 extern double p2w;
 58 extern int seg_im;
 59 
 60 extern shistogram *peaks;
 61 
 62 
 63 
 64 double nmr_hfit(int n_par, double *a, float x)
 65 {
 66    double temp = 0.0;
 67    double csf_vol,csf_grey,grey_vol,grey_white,white_vol;
 68    float nmr_wm_f_prob(double *b, float pixel_val);
 69    float nmr_gm_f_prob(double *b, float pixel_val);
 70    float nmr_csf_f_prob(double *b, float pixel_val);
 71 
 72    a[2] = a[9];
 73    if (a[1] > a[4])
 74    {
 75       a[3] = 0.0;
 76       a[4] = a[1];
 77    }
 78    if (a[4] > a[8])
 79    {
 80       a[7] = 0.0;
 81       a[8] = a[4];
 82    }
 83    temp +=  nmr_wm_f_prob(a, x) 
 84          + nmr_csf_f_prob(a, x) 
 85          + nmr_gm_f_prob(a, x);
 86    return(temp);
 87 }
 88 
 89 static double pixchisq2(int n_par, double *a)
 90 {
 91   double chisq = 0.0, temp;
 92   int     i;
 93   float x;
 94 
 95   for (i=0;i<peaks->xbins;i++)
 96     {
 97       x = peaks->xmin+(i+0.5)*peaks->xincr;
 98       temp = nmr_hfit(n_par, a , x);
 99       chisq += (peaks->array[0][i]-temp)*(peaks->array[0][i]-temp)
100         /(temp+1.0);
101     }
102   return(chisq);
103 }
104 
105 
106 void hist_fit(shistogram *hist)
107 {
108 
109     Imrect *im1;
110     Imregion *roi=NULL;
111     int n = 11;
112     double c_test1 = 0.00001;
113     double grey_vol, white_vol, csf_vol;
114     double grey_white, csf_grey;
115     double *lambda;
116 
117     if (a!=NULL) rfree(a);
118     a = (double *)ralloc(n*sizeof(double));
119     lambda = (double *)ralloc(n*sizeof(double));
120     peaks = hist;
121 
122     a[1] = csf_peak;
123     a[4] = grey_peak;
124     a[8] = white_peak;
125     a[6] = hfill1(hist,(float)((a[1]+a[4])/2.0), 0.0);
126     a[6] = sqrt(a[6]);
127     a[0] = hfill1(hist, (float)a[1], 0.0)-a[6]*a[6];
128     if (a[0]>0.0) a[0] = sqrt(a[0]);
129     else a[0] = 0.0;
130     a[3] = hfill1(hist, (float)a[4], 0.0)-a[6]*a[6];
131     if(a[3]>0.0) a[3] = sqrt(a[3]);
132     else a[3] = 0.0;
133     a[7] = hfill1(hist, (float)a[8], 0.0);
134     a[7] = sqrt(a[7]);
135     a[2] = csf_sigma;
136     a[5] = grey_sigma;
137     a[9] = white_sigma;
138     a[10] = 2.0;
139     lambda[0] = 10.0;
140     lambda[1] = 4.0;
141     lambda[2] = 1.0;
142     lambda[3] = 10.0;
143     lambda[4] = 4.0;
144     lambda[5] = 1.0;
145     lambda[6] = 10.0;
146     lambda[7] = 10.0;
147     lambda[8] = 4.0;
148     lambda[9] = 1.0;
149     lambda[10] = 5.0;
150 
151     simplexmin2(n, a, lambda, pixchisq2, c_test1, (void (*) (char *)) format);
152     hsuper(hist, n, nmr_hfit, a);
153 
154     csf_peak = a[1];
155     grey_peak = a[4];
156     white_peak = a[8];
157     csf_sigma = a[2];
158     grey_sigma = a[5];
159     white_sigma = a[9];
160     tw_fglobal_reset(pcsf_sigma);
161     tw_fglobal_reset(pgrey_sigma);
162     tw_fglobal_reset(pwhite_sigma);
163     tw_fglobal_reset(pcsf_peak);
164     tw_fglobal_reset(pgrey_peak);
165     tw_fglobal_reset(pwhite_peak);
166     csf_vol = sqrt(2.0*3.141592)*a[0]*a[0]*a[2]/hist->xincr;
167     csf_grey  = (a[4]-a[1])*a[6]*a[6]/hist->xincr;
168     grey_vol = sqrt(2.0*3.141592)*a[3]*a[3]*a[5]/hist->xincr;
169     grey_white = (a[8]-a[4])*a[10]*a[10]/hist->xincr; 
170     white_vol = sqrt(2.0*3.141592)*a[7]*a[7]*a[9]/hist->xincr;
171     format(" csf   = %3.2f + %3.2f = %3.2f \n", csf_vol,csf_grey/2.0, csf_vol+csf_grey/2.0);
172     format(" grey  = %3.2f + %3.2f + %3.2f = %3.2f \n", grey_vol, csf_grey/2.0,
173            grey_white/2.0, grey_vol + csf_grey/2.0+ grey_white/2.0);
174     format(" white = %3.2f + %3.2f = %3.2f \n", white_vol, grey_white/2.0, white_vol+grey_white/2.0);
175     rfree(lambda);
176     return;
177 }
178 
179 /*
180 routine to maximise peakiness (log entropy) of image histogram
181 using a polynomial correction to image intensity
182 void trend_proc(void)
183 {
184     Imrect *im;
185     int     type;
186     Imrect* trend_fit(Imrect *im, int m, double *am);
187     Imrect *im1;
188     void graph_hfit(Tv * tv, shistogram * hist);
189     shistogram *get_peaks_hist(void);
190     double am[4];
191  
192     if (stack_check_types(IMRECT, NULL) == false)
193     {
194         error("trend:_proc wrong type on stack", warning);
195         return;
196     }
197     im = (Imrect *) stack_pop(&type);
198     stack_push((void *)im, IMRECT,im_free);
199     am[0] = csf_peak;
200     am[1] = grey_peak;
201     am[2] = white_peak;
202     am[3] = csf_sigma;
203     im1 = trend_fit(im, 4, (double *)am);
204     stack_push((void *) im1, IMRECT, im_free);
205     imcalc_draw(imcalc_tv_get());
206     imcalc_draw(imcal2_tv_get());
207     image_choice_reset();
208     graph_hfit(imcalc_graph_tv_get(), get_peaks_hist());
209 }
210 */
211 
212 
213 void     nmr_fit_proc(void)
214 {
215     Imrect *im;
216     int     type;
217     char   temp[1024];
218     Imregion *roi;
219     Imrect *im_hist(shistogram **imhist, double k, double range,
220                     Imrect *im, Imregion *roi);
221     void hist_fit(shistogram *hist);
222     void graph_hfit(Tv * tv, shistogram * hist);
223     double range, k;
224     Tv *imcalc_graph_tv_get();
225 
226     if (stack_check_types(IMRECT, NULL) == false)
227     {
228         error("imcalc_roi: wrong type on stack", warning);
229         return;
230     }
231     im = (Imrect *) stack_inspect(&type);
232 
233     if (peaks != NULL)
234       hfree(peaks);
235     peaks = NULL;
236 
237 /*
238     k = (csf_peak-3.0*csf_sigma + white_peak+3.0*white_sigma)/2.0;
239     range = (white_peak+3.0*white_sigma - csf_peak+3.0*csf_sigma)/2.0;
240 */
241     k = (csf_min + white_max) /2.0;
242     range = (white_max - csf_min ) /2.0;
243 
244     roi = tv_get_im_roi(imcalc_tv_get()); 
245     im_free(im_hist(&peaks, k, range, im, roi));
246     sprintf(temp, "overflows %f underflows %f entries %f mean %f \n",
247             peaks->over,peaks->under,(float)peaks->entries,(float)peaks->mean);
248     format(temp);
249     hist_fit(peaks);
250     graph_hfit(imcalc_graph_tv_get(), peaks);
251 }
252 
253 void gaus_and_tri_proc()
254 {
255     void gaus_and_blur_tri_func(double *newpar, int select);
256 
257     if (a!=NULL) gaus_and_blur_tri_func(a,seg_im);
258 }
259 
260 
261 void    rusenik(void)
262 {
263     Imrect *im,*im2;
264     int     type;
265     Imregion *roi;
266     Complex *row1;
267     float *row2;
268     int     lx, ux, ly, uy;
269     int     i, j;
270     int width, height;
271     double  s1,s2;
272     double  c=0,w=0,g=0;
273     double c_tot=0.0,w_tot=0.0,g_tot=0.0;
274 
275     if (stack_check_types(IMRECT, NULL) == false)
276     {
277         error("imcalc_roi: wrong type on stack", warning);
278         return;
279     }
280     im = (Imrect *) stack_inspect(&type);
281     if (im->vtype!=complex_v) return;
282 
283     roi = im->region;
284     if (roi == NULL)
285         return;
286     width = im->width;
287     height = im->height;
288 
289     lx = roi->lx;
290     ux = roi->ux;
291     ly = roi->ly;
292     uy = roi->uy;
293 
294     im2 = im_alloc(height, width, roi, float_v);
295     row1 = zvector_alloc(lx, ux);
296     row2 = fvector_alloc(lx, ux);
297 
298 
299     for (i = ly; i < uy; ++i)
300     {
301         im_get_rowz(row1, im, i, lx, ux);
302         for (j = lx; j < ux; ++j)
303         {
304             s1 = row1[j].x;
305             s2 = row1[j].y;
306             if (s1!=0.0 && s2 !=0.0)
307             {
308                c = (s1*(p2w-p2g) - s2*(p1w-p1g) - (p1g*p2w - p2g*p1w))/
309                  ( (p1c-p1g)*(p2w-p2g) - (p2c-p2g)*(p1w-p1g) ); 
310 
311                w = (s1*(p2c-p2g) - s2*(p1c-p1g) - (p1g*p2c - p2g*p1c))/
312                  ( (p2c-p2g)*(p1w-p1g) - (p1c-p1g)*(p2w-p2g) );
313                g = 1.0 - c - w;
314                if (seg_im==0) row2[j] = c;
315                if (seg_im==1) row2[j] = g;
316                if (seg_im==2) row2[j] = w;
317             }
318             else row2[j] = 0.0;
319             c_tot += c;
320             w_tot += w;
321             g_tot += g;
322         }
323         im_put_rowf(row2, im2, i, lx, ux);
324     }
325     format("csf = %3.2f \n",c_tot);
326     format("grey = %3.2f \n",g_tot);
327     format("white = %3.2f \n",w_tot);
328 
329     zvector_free( row1, lx);
330     fvector_free( row2, lx);
331     stack_push(im2,IMRECT, im_free);
332     imcalc_draw(imcalc_tv_get());
333     imcalc_draw(imcal2_tv_get());
334     image_choice_reset();
335 }
336 

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