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

Linux Cross Reference
Tina6/tina-libs/tina/image/imgEM_grad.c

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

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

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