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

Linux Cross Reference
Tina5/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 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 

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