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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.