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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedSeg_scale.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    :  $Source: /home/tina/cvs/tina-tools/tinatool/tlmedical/tlmedSeg_scale.c,v $
 37  * Date    :  $Date: 2004/02/13 12:41:23 $
 38  * Version :  $Revision: 1.4 $
 39  * CVS Id  :  $Id: tlmedSeg_scale.c,v 1.4 2004/02/13 12:41:23 neil Exp $
 40  *
 41  * Notes :
 42  *
 43 */
 44 
 45 #include "tlmedSeg_scale.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 static double *a = NULL;
 62 static shistogram *peaks;
 63 static Mixmodel *model_1d =NULL;
 64 
 65 
 66 static double nmr_hfit(int n_par, double *a, float x)
 67 {
 68    double temp = 0.0;
 69    int i,j,k;
 70 
 71    k = 0;
 72    for (i=0;i<model_1d->nmix;i++)
 73    {
 74        for (j=i;j<model_1d->nmix;j++)
 75        {
 76            if (model_1d->par_dens[i][j]>0)
 77            {
 78                k++;
 79                model_1d->par_dens[i][j] = fabs(a[k]);
 80                model_1d->par_dens[j][i] = fabs(a[k]);
 81            }
 82        }
 83    }
 84    for (i=0;i<model_1d->nmix;i++)
 85        model_1d->normal[i] = model_1d->par_dens[i][i]; 
 86 
 87    temp = em_graph_data(model_1d, a[0]*x);
 88    return(temp*peaks->xincr);
 89 }
 90 
 91 static double pixchisq2(int n_par, double *a)
 92 {
 93   double chisq = 0.0, temp;
 94   int     i;
 95   float x;
 96 
 97   for (i=0;i<peaks->xbins;i++)
 98     {
 99       x = peaks->xmin+(i+0.5)*peaks->xincr;
100       temp = nmr_hfit(n_par, a , x);
101       chisq += (peaks->array[0][i]-temp)*(peaks->array[0][i]-temp)
102         /(temp+1.0);
103     }
104   return(chisq);
105 }
106 
107 static void hist_fit(shistogram *hist, int imno, Mixmodel *model)
108 {
109     int n = 1;
110     double c_test1 = 0.00001;
111     double *lambda;
112     int i,j,k;
113     
114     model_1d = get_submodel(model, imno , imno );
115     if (a!=NULL) rfree(a);
116     for (i=0;i<model->nmix;i++)
117     {
118         for (j=i;j<model->nmix;j++)
119         {
120             if (model->par_dens[i][j]>0)
121             {
122                 n++;
123             }
124         }
125     }
126 
127     a = (double *)ralloc(n*sizeof(double));
128     lambda = (double *)ralloc(n*sizeof(double));
129     
130     /* Image Scale */  
131     a[0] = 1.0; 
132     lambda[0] = 0.02;
133 
134     /* normalisations */
135     k = 0;
136     for (i=0;i<model->nmix;i++)
137     {
138         for (j=i;j<model->nmix;j++)
139         {
140             if (model->par_dens[i][j]>0)
141             {
142                 k++;
143                 a[k] = model->par_dens[i][j];
144                 lambda[k] = sqrt(a[k]);
145             }
146         }
147     }
148 
149     for(i = 0;i < n;i++)
150        format("a%i = %f  ", i,a[i]);
151 
152     simplexmin2(n, a, lambda, pixchisq2, NULL, c_test1, (void (*) (char *)) format);
153     hsuper(hist, n, nmr_hfit, a);
154     for(i = 0;i < n;i++)
155        format("a%i = %f  ", i,a[i]);
156 
157     for (i=0;i<model->nmix;i++) model->vectors[i][imno] /= a[0]; 
158 
159     for (i=0;i<model->nmix;i++) model->cov[i][imno][imno] /= 1.0/(a[0]*a[0]);
160 
161     k = 0;
162     for (i=0;i<model->nmix;i++)
163     {
164         for (j=i;j<model->nmix;j++)
165         {
166             if (model->par_dens[i][j]>0)
167             {
168                 k++;
169                 model->par_dens[i][j] = fabs(a[k]);
170                 model->par_dens[j][i] = fabs(a[k]);
171             }
172         }
173     }
174     for (i=0;i<model->nmix;i++)
175          model->normal[i] = model->par_dens[i][i]; 
176 
177     rfree(lambda);
178     return;
179 }
180 
181 shistogram    *nmr_scaled_fit(Imregion *roi, Mixmodel *m)
182 {
183     Imrect *im;
184     int im_no,xmax,ymax,xmin,ymin;
185     char temp[1024];
186     double max_g,min_g;   
187     double range, k;
188     List *lptr;
189     Sequence *seq = seq_get_current();
190 
191     if(seq == NULL)
192     {
193        error("There is no sequence! \n", warning);
194        return(NULL);
195     }
196     if(m == NULL)
197     {
198        error("Model parameters not defined!\n", warning);
199        return(NULL);
200     }
201 
202     lptr = get_seq_current_el(seq);
203     im = (Imrect *)(lptr->to);
204     im_no = get_seq_current_frame(seq);
205 
206     /* Check if number of images in sequence are more than number of parameters specified */
207 /*
208     if((seq->seq_end - seq->seq_start) > m->ndim)
209        im_no = 0;
210 */
211 
212     max_g = imf_max(im, &ymax,&xmax);
213     min_g = imf_min(im, &ymin,&xmin);
214 
215     if (peaks != NULL)
216       hfree(peaks);
217     peaks = NULL;
218 
219     k = (min_g + max_g)/2;
220     range = fabs(max_g-min_g)/2; 
221   
222     im_free(im_hist(&peaks, k, range, im, roi));
223     sprintf(temp, "overflows %f underflows %f entries %f mean %f \n",
224             peaks->over,peaks->under,(float)peaks->entries,(float)peaks->mean);
225     format(temp);
226     hist_fit(peaks,im_no,m);
227     return(peaks);
228 }
229 

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