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 Lesser 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 Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser 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 **********
20 *
21 * Program : TINA
22 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgEM_grad.c,v $
23 * Date : $Date: 2005/02/07 23:27:41 $
24 * Version : $Revision $
25 * CVS Id : $Id: imgEM_grad.c,v 1.4 2005/02/07 23:27:41 paul Exp $
26 *
27 * Notes :
28 *
29 *********
30 */
31
32 #include "imgEM_grad.h"
33
34 #if HAVE_CONFIG_H
35 # include <config.h>
36 #endif
37
38 #include <stdio.h>
39 #include <math.h>
40 #include <tina/sys/sysDef.h>
41 #include <tina/sys/sysPro.h>
42 #include <tina/math/mathDef.h>
43 #include <tina/math/mathPro.h>
44 #include <tina/image/imgDef.h>
45 #include <tina/image/imgPro.h>
46
47 #define MAXPLOT 255
48 #define OFFSET 0.35
49 #define GAMMA 1
50
51 extern double **get_kgrad_arr(void);
52
53 /*-------------------------------------------------*/
54 /* Range of values for scaling original grey level */
55 /* image and associated gradient image */
56 /*-------------------------------------------------*/
57 static double scale_grey = 255.0;
58 static double scale_slope = 255.0;
59 static double low_slope = 0.0;
60 static double low_grey = 0.0;
61 static double *noise_arr = NULL;
62
63 Imrect *mD_slope_image(Sequence *seq)
64 {
65 Imrect *im_tot=NULL, *im_tot1=NULL, *im_seq=NULL, *im=NULL;
66 Imrect *im1=NULL, *im2=NULL, *im3=NULL, *im_norm=NULL;
67 Imregion *roi = NULL;
68 int lx, ux, ly, uy;
69 List *im_ptr;
70 double noise_k,sigma_x, sigma_y, diffxy;
71 int count = 0, k=0;
72
73 if (seq == NULL)
74 {
75 error("Sequence is missing!", warning);
76 return(NULL);
77 }
78 if ((im_ptr = get_seq_start_el(seq)) == NULL)
79 return(NULL);
80
81 while(im_ptr->next != NULL)
82 {
83 im_ptr = im_ptr->next;
84 count++;
85 }
86 if (noise_arr !=NULL) dvector_free(noise_arr,0);
87 noise_arr = dvector_alloc(0, count+1);
88
89 im_ptr = get_seq_start_el(seq);
90 im_seq = (Imrect *)im_ptr->to;
91 roi = im_seq->region;
92 if (roi == NULL)
93 return(NULL);
94 lx = roi->lx;
95 ux = roi->ux;
96 ly = roi->ly;
97 uy = roi->uy;
98 im_tot = (void *) im_alloc(uy-ly,ux-lx,roi,float_v);
99
100 while(im_ptr->to != NULL)
101 {
102 im_seq = (Imrect *)im_ptr->to;
103 im1 = imf_diffx(im_seq);
104 im2 = im_sqr(im1);
105 im_free(im1);
106
107 im1 = imf_diffy(im_seq);
108 im3 = im_sqr(im1);
109 im_free(im1);
110
111 im = im_sum(im2,im3);
112 /* Noise calculation as in "noise" Imcacl button */
113 sigma_x = imf_diffx_noise(im_seq,roi);
114 sigma_y = imf_diffy_noise(im_seq,roi);
115 diffxy = fabs(sigma_x-sigma_y);
116 noise_k = (sigma_x + sigma_y)/2.0;
117 noise_arr[k] = noise_k;
118 k++;
119 /* Normalize gradient image to noise in original image*/
120 im_norm = im_times(1.0/(noise_k*noise_k),im);
121 im_free(im);
122 im_tot1 = im_sum(im_tot,im_norm);
123 im_free(im_tot);
124 im_free(im_norm);
125 im_tot = im_tot1;
126 im_ptr = im_ptr->next;
127 if (im_ptr == NULL) break;
128 }
129 im_free(im2);
130 im_free(im3);
131 im = im_sqrt(im_tot);
132 if (count > 1)
133 {
134 im2 = imf_add( OFFSET * (sqrt(count)-count),im );
135 im_free(im);
136 im = im2;
137 }
138 im_free(im_tot);
139 return(im);
140 }
141
142
143 /*************************************************************************/
144 /* Generate scatter plot Gradient vs Grey Level normalized to noise */
145 /*************************************************************************/
146
147 Imrect *mD_scat_grad(Imrect *im, Imrect *im_slope)
148 {
149 Imrect *im_scat = NULL, *im_scale = NULL, *im_slope_scale = NULL;
150 Imrect *im_norm;
151 Imregion *roi_scat = NULL, *roi=NULL;
152 int lx, ux, ly, uy;
153 int i,j,x,y;
154 double pixval, range_x, range_y;
155 float *rowx, *rowy;
156 double noise_k,sigma_x, sigma_y, diffxy;
157
158 if (im == NULL) return(NULL);
159 if (im_slope == NULL) return(NULL);
160
161 if (roi == NULL ) roi = im->region;
162 lx = roi->lx;
163 ux = roi->ux;
164 ly = roi->ly;
165 uy = roi->uy;
166 format(" lx = %i ux = %i ly = %i uy = %i\n",lx,ux,ly,uy);
167
168 /* Noise calculation as in "noise" Imcalc button */
169 sigma_x = imf_diffx_noise(im,roi);
170 sigma_y = imf_diffy_noise(im,roi);
171 diffxy = fabs(sigma_x-sigma_y);
172 noise_k = (sigma_x + sigma_y)/2.0;
173
174 /* Normalize original image to noise in original image*/
175 im_norm = im_times(1.0/noise_k,im);
176
177 im_scale = imf_scale(im_norm, low_grey , scale_grey);
178 im_slope_scale = imf_scale(im_slope, low_slope, scale_slope);
179
180 if (im_scale == NULL) return (NULL);
181 if (im_slope_scale == NULL) return(NULL);
182
183 range_x =(MAXPLOT+1.0)/scale_grey;
184 range_y =(MAXPLOT+1.0)/scale_slope;
185
186 format(" range_x = %f range_y = %f\n",range_x,range_y);
187
188 /*im_scat = im_alloc(im->height, im->width, NULL, float_v);*/
189 /*im_scat = im_alloc(uy-ly,ux-lx,roi,float_v);*/
190 im_scat = im_alloc(256,256,NULL,float_v);
191 roi_scat = im_scat->region;
192
193 rowx = fvector_alloc(lx,ux);
194 rowy = fvector_alloc(lx,ux);
195
196 for( i = ly+1; i < uy-1; ++i)
197 {
198 im_get_rowf(rowx,im_scale,i,lx,ux);
199 im_get_rowf(rowy,im_slope_scale,i,lx,ux);
200 for (j = lx+1; j < ux-1; ++j)
201 {
202 x = (int) rowx[j];
203 y = (int) rowy[j];
204 if( (int)(range_x*x) > lx && (int)(range_x*x) < ux &&
205 (int)(range_y*y) > ly && (int)(range_y*y) <uy )
206 {
207 pixval = im_get_pix(im_scat, MAXPLOT - (int)(range_y* y),(int)(range_x*x));
208 pixval += 1.0;
209 im_put_pix(pixval, im_scat,MAXPLOT - (int)(range_y* y),(int) (range_x*x));
210 }
211 else
212 im_put_pix(0.0, im_scat,MAXPLOT - (int)(range_y* y),(int) (range_x*x));
213 }
214 }
215 fvector_free(rowx, lx);
216 fvector_free(rowy, lx);
217 im_free(im_scale);
218 im_free(im_slope_scale);
219 im_free(im_norm);
220 return(im_scat);
221 }
222
223
224
225 double mD_dist_calc(Mixmodel *m, int blob1, int blob2)
226 {
227 int i,n;
228 double *m1, *m2;
229 double d = 0.0, r=0.0;
230
231 if (m == NULL)
232 {
233 error("Input file is missing!", warning);
234 return(0);
235 }
236 if(noise_arr == NULL)
237 {
238 error("Noise needs to be calculated!",warning);
239 return(0);
240 }
241 if (blob1 == blob2) return(0.0);
242 n = m->ndim;
243 m1 = m->vectors[blob1];
244 m2 = m->vectors[blob2];
245 for(i = 0; i < n; i++)
246 d += (m1[i]-m2[i])*(m1[i]-m2[i])/(noise_arr[i]*noise_arr[i]);
247 r = sqrt(d);
248 return(r);
249 }
250
251 double mD_grad_func(double f, double d, double y, int loop, int loop1, double **sk, double s0, double *sf)
252 {
253 double Pgrad=0.0,w,s,news0;
254 /* x - grey level vaulues */
255 /* y - gradient */
256
257 if (sk[loop][loop1]==0.0) return(0.0);
258 if (y < 0.001) y = 0.001;
259
260 if(loop == loop1)
261 *sf = sk[loop][loop1]*s0;
262 else
263 {
264 /* Calculate weighting factor for distance */
265 w = 1.0 - 4.0*(f-0.5)*(f-0.5);
266 /* Calculate sigma(f) */
267 news0 = s0*(f*sk[loop][loop] + (1.0-f)*sk[loop1][loop1])/sk[loop][loop1];
268
269 *sf = sk[loop][loop1]*sqrt(news0*news0 + w*d*d/4.0);
270 }
271 /* Calculate gradiant distribution value Pgrad for gradient value corresponding to current gray level value */
272 s = *sf * *sf;
273 if (GAMMA==1) Pgrad = y*exp(-0.5*(y*y/(s)))/s;
274 else Pgrad = y*y*exp(-0.5*(y*y/(s)))/(s* *sf);
275
276 return(Pgrad);
277 }
278
279
280 Imrect *mD_grad_model_norm(Imrect *im,Imrect *im_slope, Mixmodel *m, Mixmodel *m_sub)
281 {
282 double s0=1.0,Pgrad=0.0,sf=0.0;
283 int n;
284 int lx, ux, ly, uy;
285 int i,j, loop, loop1;
286 Imrect *im_model;
287 Imregion *roi;
288 float *row;
289 float x[1];
290 double y;
291 double high=256.0, low1=0.0, low2=0.0, scale1, scale2;
292 double frac1, frac2;
293 float min1, max1, min2, max2;
294 double **Sk,a,d;
295
296
297 if (im == NULL)
298 {
299 error("Sequence image is missing!", warning);
300 return(NULL);
301 }
302 if (im_slope == NULL)
303 {
304 error("Gradient image is missing!", warning);
305 return(NULL);
306 }
307 roi = im_slope->region;
308 if (roi == NULL)
309 return(NULL);
310 if (m == NULL)
311 {
312 error("Input file is missing!", warning);
313 return(NULL);
314 }
315 lx = roi->lx;
316 ux = roi->ux;
317 ly = roi->ly;
318 uy = roi->uy;
319 im_model = (void *) im_alloc(uy-ly,ux-lx,roi,float_v);
320 row = fvector_alloc(lx, ux);
321
322 imf_minmax(im, &min1, &max1);
323 if (max1 == min1)
324 max1 = (double) (min1 + 1.0);
325 scale1 = (double) fabs((high - low1) / (max1 - min1));
326 low1 -= (double)(scale1 * min1);
327
328 imf_minmax(im_slope, &min2, &max2);
329 if (max2 == min2)
330 max2 = (double) (min2 + 1.0);
331 scale2 = (double)fabs((high - low2) / (max2 - min2));
332 low2 -= (double)(scale2 * min2);
333
334 /* Calculate sigma0 */
335 n = m->ndim;
336 s0 = (double)sqrt(n);
337 /* Get the values of scaling factors Sk */
338 Sk = m->gradscale;
339 if(Sk == NULL)
340 {
341 error("Calculate gradient scaling parameters <k grad calc>", warning);
342 return(NULL);
343 }
344 /* Calculate normalisation factors for gradient distribution */
345
346 for (i = 0; i < high; ++i) /* gradient loop */
347 {
348 for (j = 0; j < high; ++j) /* grey level loop */
349 {
350 x[0] = (float)(j + 0.5 -low1)/scale1;
351 y = (double)(i + 0.5 -low2)/scale2;
352 /* Loop over tissue models */
353 for(loop = 0; loop < m->nmix; loop++)
354 {
355 for(loop1 = 0; loop1 < m->nmix; loop1++)
356 {
357 if (loop == loop1)
358 {
359 Pgrad += m->normal[loop]*pure_density(m_sub,(float*)&x[0],loop)*
360 mD_grad_func(1.0,0.0,y,loop,loop1,Sk,s0,&sf);
361 }
362 else
363 {
364 a = part_fraction(m_sub,(float*)&x[0],loop,loop1);
365 d = mD_dist_calc(m,loop,loop1);
366 frac1 = part_density(m_sub,(float*)&x[0],loop,loop1,a);
367 frac2 = part_density(m_sub,(float*)&x[0],loop1,loop,(1.0-a));
368 if (frac1+frac2 != 0.0) a = frac1/(frac1+frac2);
369 Pgrad += m->par_dens[loop][loop1]*frac1*
370 mD_grad_func(a,d,y,loop,loop1,Sk,s0,&sf);
371 if(Pgrad >=0.0 || Pgrad <=0.0)
372 {
373 }
374 else
375 {
376 printf("nan\n");
377 }
378
379 }
380 }
381 }
382 row[j] = (float) Pgrad;
383 Pgrad =0.0;
384 }
385 im_put_rowf(row,im_model,MAXPLOT-i, lx, ux);
386 }
387 fvector_free(row,lx);
388
389 return(im_model);
390 }
391
392 double grad_hfit(int n, double *a, float x)
393 {
394 double temp;
395
396 if (GAMMA ==1)
397 temp = 0.177*a[0]*x*exp(-x*x/2.0);
398 else
399 temp = 0.177*a[0]*x*x*exp(-x*x/2.0);
400 return (temp);
401 }
402
403
404 double **mD_grad_kcalc(Mixmodel *m, Imrect *im, Imrect *im_slope,void ***pim_array, void ***pim_array_grad)
405 {
406 shistogram *imhist1,*imhist2;
407 double *am, *bm, f, d;
408 void ***imptrs = NULL;
409 int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
410 int lx, ux, ly, uy;
411 int i,j,k,loop,loop1;
412 float *row, *data;
413 Imregion *roi;
414 List *lptr;
415 Sequence *seq = seq_get_current();
416 double **Sk, y, grad_mean;
417 double s0=1.0,sf=0.0;
418 Imrect *im_tot1=NULL, *im1=NULL, *im2=NULL;
419 double tot, this_pix1, this_pix2, **M1, **M2;
420 float grad_min,grad_max;
421
422 if (m == NULL)
423 {
424 error("Input file is missing!", warning);
425 return(NULL);
426 }
427 if (seq == NULL)
428 {
429 error("Sequence is missing!", warning);
430 return(NULL);
431 }
432 if (pim_array_grad == NULL)
433 {
434 error("Using grey levels only for probability calculation.", warning);
435 }
436
437 imhist1 = hbook1(0,"gradient pure\0",0.0, 5.0, 25);
438 imhist2 = hbook1(0,"gradient mix\0",0.0, 5.0, 25);
439
440
441 lptr = get_seq_current_el(seq);
442 im = lptr->to;
443 roi = im_slope->region;
444
445 lx = roi->lx;
446 ux = roi->ux;
447 ly = roi->ly;
448 uy = roi->uy;
449
450 seq_slice_init(seq);
451 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
452 &imptrux);
453 data = (float *)fvector_alloc(imptrlz, imptruz);
454
455 M1 = darray_alloc(0, 0, m->nmix, m->nmix);
456 M2 = darray_alloc(0, 0, m->nmix, m->nmix);
457 imf_minmax(im_slope, &grad_min, &grad_max);
458 Sk = m->gradscale;
459 if(Sk == NULL)
460 {
461 Sk = darray_alloc(0,0,m->nmix,m->nmix);
462 for(loop = 0; loop < m->nmix; loop++)
463 for(loop1 = 0; loop1 < m->nmix; loop1++)
464 {
465 if (loop==loop1) Sk[loop][loop1] = 1.0;
466 else Sk[loop][loop1] = 0.5;
467 }
468 }
469 s0 = (double)sqrt(m->ndim);
470 /* Loop over all tissues */
471 for(loop = 0; loop < m->nmix; loop++)
472 {
473 for(loop1 = 0; loop1 < m->nmix; loop1++)
474 {
475 if (m->par_dens[loop][loop1] == 0.0) continue;
476 im1 = (Imrect *) pim_array[loop][loop1];
477 im2 = (Imrect *) pim_array[loop1][loop];
478 if(im_tot1!=NULL)
479 {
480 im_free(im_tot1); /* Fixes major memory leak: im_tot1 is allocated inside this loop, */
481 im_tot1=NULL; /* but was being freed outside it PAB 7/2/2005 */
482 }
483 im_tot1 = prob_im_tot(m,pim_array,roi);
484 M1[loop][loop1] = 0.0;
485 M2[loop][loop1] = 0.0;
486 d = mD_dist_calc(m,loop,loop1);
487 /* Loop over all data */
488 for (i = ly; i < uy; ++i)
489 {
490 for (j = lx; j < ux; ++j)
491 {
492 for (k = imptrlz; k < imptruz; k++)
493 {
494 row = imptrs[k][i];
495 data[k] = row[j];
496 }
497 /* Calculate mean of the data */
498 y = im_get_pixf(im_slope, i, j);
499 this_pix1 = im_get_pixf(im1,i,j);
500 this_pix2 = im_get_pixf(im2,i,j);
501 tot = im_get_pixf(im_tot1,i,j);
502 if (this_pix1+this_pix2!=0)
503 f = this_pix1/(this_pix1+this_pix2);
504 else
505 f = 0.0;
506 mD_grad_func(f,d,y,loop,loop1,Sk,s0,&sf);
507
508 if (GAMMA==1) grad_mean = sqrt(PI/2)*sf;
509 else grad_mean = sqrt(8.0/PI)*sf;
510
511 if (tot > 0.0 && y/sf <= 4.0)
512 {
513 M1[loop][loop1] += y*this_pix1/tot;
514 M2[loop][loop1] += grad_mean*this_pix1/tot;
515 if (loop == loop1)
516 {
517 hfill1(imhist1,y/sf,this_pix1/tot);
518 }
519 else
520 {
521 hfill1(imhist2,y/sf,this_pix1/tot);
522 }
523 }
524 }
525 }
526 }
527 }
528 for (loop = 0; loop < m->nmix; loop++)
529 for (loop1 = 0; loop1 < m->nmix; loop1++)
530 {
531 if ( (M2[loop][loop1]) > 0.0 )
532 Sk[loop][loop1] *= (M1[loop][loop1]+M1[loop1][loop])
533 /(M2[loop][loop1]+M2[loop1][loop]);
534 format("i = %i, j = %i, Sk = %f \n",loop,loop1,Sk[loop][loop1]);
535 }
536
537 im_free(im_tot1);
538 im_tot1=NULL;
539 fvector_free(data,imptrlz);
540 darray_free(M1, 0, 0, m->nmix, m->nmix);
541 darray_free(M2, 0, 0, m->nmix, m->nmix);
542 am = dvector_alloc(0,4);
543 am[0] = imhist1->contents;
544 bm = dvector_alloc(0,4);
545 bm[0] = imhist2->contents;
546 hsuper(imhist1, 4, grad_hfit, am);
547 hsuper(imhist2, 4, grad_hfit, bm);
548 hprint(stdout,imhist1);
549 hprint(stdout,imhist2);
550 return(Sk);
551 }
552
553
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.