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/medical/medNorm_seq.c,v $
23 * Date : $Date: 2004/12/06 22:05:08 $
24 * Version : $Revision: 1.5 $
25 * CVS Id : $Id: medNorm_seq.c,v 1.5 2004/12/06 22:05:08 paul Exp $
26 *
27 * Notes :
28 *
29 * Gio Buonaccorsi removed the #if 0 and altered seq_norm to return a list
30 * of imrects. Purpose was to improve lib / tool separation.
31 * Stack push of imrects has been moved to the calling function, batch_norm_proc in
32 * tlmedCoreg_tool.c.
33 *
34 *********
35 */
36
37 #include "medNorm_seq.h"
38
39 #if HAVE_CONFIG_H
40 # include <config.h>
41 #endif
42
43 #include <stdio.h>
44 #include <math.h>
45 #include <tina/sys/sysDef.h>
46 #include <tina/sys/sysPro.h>
47 #include <tina/math/mathDef.h>
48 #include <tina/math/mathPro.h>
49 #include <tina/image/imgDef.h>
50 #include <tina/image/imgPro.h>
51 #include <tina/medical/medNorm_base.h>
52
53 /* A reasonable attempt at 3D coil correction based upon the 2D algorithm */
54 /* assumes approx equal errors on log xy correction field */
55 /* Gaussian smoothing in z approximated by repeated convolution with */
56 /* a double sided exponential filter */
57 /* NAT 5/5/2002 */
58
59 static void vol_smooth(Imrect **ims, double sigma, int data_start, int data_end)
60 /* double sided exponential smoothing as in xy correction algorithm */
61 /* with edge correction field fixed */
62 {
63 int i;
64 double a, b = exp(-1.0 / sigma), ab;
65 Imrect *im1;
66 Imrect **impl, **immi;
67
68 a = (1 - b)/( 1 + b );
69 ab = a*b;
70
71 impl = (Imrect **)pvector_alloc(data_start,data_end+1);
72 immi = (Imrect **)pvector_alloc(data_start,data_end+1);
73
74 for (i=data_start; i<=data_end; i++)
75 {
76 if ( i == data_start)
77 im1 = imf_wsum(a,b,ims[i],ims[i-1]);
78 else
79 im1 = imf_wsum(a,b,ims[i],impl[i-1]);
80 impl[i] = im1;
81 }
82 for (i=data_end; i>=data_start; i--)
83 {
84 if (i == data_end)
85 im1 = imf_wsum(ab,0.0,ims[i+1],ims[i+1]);
86 else
87 im1 = imf_wsum(ab,b,ims[i+1],immi[i+1]);
88 immi[i] = im1;
89 }
90 for (i=data_start; i<=data_end; i++)
91 {
92 im1 = imf_sum(impl[i],immi[i]);
93 im_free(ims[i]);
94 im_free(impl[i]);
95 im_free(immi[i]);
96 ims[i] = im1;
97 }
98 pvector_free(impl,data_start);
99 pvector_free(immi,data_start);
100 }
101
102 /* Static but not used: comment out for now PAB */
103 /* double sided exponential smoothing as in xy correction algorithm */
104 /* with edge correction field fixed */
105 /*
106 static void vol_smooth2(Imrect **ims, double sigma, int data_start, int data_end)
107 {
108 int i;
109 Imrect *im1;
110 Imrect **impl, **immi;
111
112 impl = (Imrect **)pvector_alloc(data_start,data_end+1);
113 immi = (Imrect **)pvector_alloc(data_start,data_end+1);
114
115 for (i=data_start; i<=data_end; i++)
116 {
117 im1 = imf_wsum(0.5,0.5,ims[i],ims[i-1]);
118 impl[i] = im1;
119 }
120 for (i=data_end; i>=data_start; i--)
121 {
122 im1 = imf_wsum(0.5,0.5,ims[i],ims[i+1]);
123 immi[i] = im1;
124 }
125 for (i=data_start; i<=data_end; i++)
126 {
127 im1 = imf_wsum(0.5,0.5,impl[i],immi[i]);
128 im_free(ims[i]);
129 im_free(impl[i]);
130 im_free(immi[i]);
131 ims[i] = im1;
132 }
133 pvector_free(impl,data_start);
134 pvector_free(immi,data_start);
135 }
136 */
137
138 static void z_smooth(double *factor, double sigma, int data_start, int data_end)
139 {
140 double *factorf, *factorb;
141 double a, b = exp(-1.0 / sigma), ab;
142 int i;
143
144 a = (1 - b)/( 1 + b );
145 ab = a*b;
146
147 factorf = (double *)dvector_alloc(data_start-1,data_end+2);
148 factorb = (double *)dvector_alloc(data_start-1,data_end+2);
149 for (i=data_start; i<=data_end; i++)
150 {
151 if (i == data_start)
152 factorf[i] = a*factor[i]+b*factor[i-1];
153 else
154 factorf[i] = a*factor[i]+b*factorf[i-1];
155 }
156 for (i=data_end; i>=data_start; i--)
157 {
158 if (i == data_end)
159 factorb[i] = ab*factor[i+1]+b*factor[i+1];
160 else
161 factorb[i] = ab*factor[i+1]+b*factorb[i+1];
162 }
163 for (i=data_start; i<=data_end; i++)
164 {
165 factor[i] = factorf[i]+factorb[i];
166 }
167 dvector_free(factorf, data_start-1);
168 dvector_free(factorb, data_start-1);
169 }
170
171 /*
172 * Added mod for tool / lib separation - GAB 18 Feb 2003
173 *
174 * Altered to return image list (for stack push in caller).
175 */
176 List *seq_norm(Sequence *data, double sigma, double thresh)
177 {
178 int i, seq_start, seq_end;
179 List *lptr, *imlistptr = NULL;
180 Imrect *im, *im1, *im2, *cut, *image1, *image2, **ims;
181 double noise = 0.0, *factor;
182
183 if (data==NULL) return NULL;
184 seq_start = data->offset;
185 seq_end = get_end_frame(data);
186
187
188 ims = (Imrect **)pvector_alloc(seq_start-1,seq_end+2);
189
190 for (i=seq_start; i<=seq_end; i++)
191 {
192 set_seq_current_frame(data, i);
193 lptr = get_seq_current_el(data);
194 if ((lptr == NULL) || (im = lptr->to) == NULL) break;
195 format("xy scale image %d \n", i);
196 if (noise == 0.0)
197 {
198 noise = imf_diffx_noise(im, NULL);
199 noise += imf_diffy_noise(im, NULL);
200 noise *= 0.5;
201 format("estimated noise %f \n",noise);
202 }
203 /* calculate logarithm of the xy correction field */
204 ims[i] = xy_norm(im, noise, sigma, thresh);
205 }
206 /* pad out start and end for smoothing */
207 ims[i] = im_copy(ims[i-1]);
208 ims[seq_start-1] = im_copy(ims[seq_start]);
209
210 /* exponential smoothing assuming equal errors on log correction */
211 vol_smooth(ims, 1.0, seq_start, seq_end);
212 vol_smooth(ims, 1.0, seq_start, seq_end);
213
214 /* convert to multiplicative xy field */
215 for (i=seq_start; i<=seq_end; i++)
216 {
217 im1 = imf_exp(ims[i]);
218 im_free(ims[i]);
219 ims[i] = im1;
220 }
221
222 set_seq_current_frame(data, seq_start);
223 lptr = get_seq_current_el(data);
224 im1 = im_prod(lptr->to,ims[seq_start]);
225 factor = (double *)dvector_alloc(seq_start-1,seq_end+2);
226 factor[seq_start] = 0.0;
227
228 for (i=seq_start+1; i<=seq_end; i++)
229 {
230 set_seq_current_frame(data, i);
231 lptr = get_seq_current_el(data);
232 if ((lptr == NULL) || (im = lptr->to) == NULL) break;
233 im2 = im_prod(im,ims[i]);
234
235 edge_mask(noise,im1,im2,&cut,&image1,&image2);
236 factor[i] = im_fraction(cut,image1,image2,31,noise);
237 format("z scale factor %d = %f \n", i, factor[i]);
238 im_free(im1);
239 im_free(cut);
240 im_free(image1);
241 im_free(image2);
242 im1 = im2;
243 }
244 /* pad out start and end for smoothing */
245 factor[seq_start-1] = 0.0;
246 factor[seq_end+1] = factor[seq_end];
247 /* exponential smoothing preserving integral product across z map*/
248 z_smooth(factor, 1.5, seq_start, seq_end);
249 z_smooth(factor, 1.5, seq_start, seq_end);
250
251 factor[seq_start-1] = 1.0;
252 for (i=seq_start; i<=seq_end; i++)
253 {
254 factor[i] = exp(factor[i]);
255 if (i > seq_start) factor[i] *= factor[i-1];
256 format("z scale factor %d = %f \n", i, factor[i]);
257
258 imf_times_inplace(factor[i],ims[i]);
259
260 /*
261 * Attempt at add to list approach - GAB 18 Feb 2003
262 */
263 imlistptr = ref_addtostart(imlistptr, ims[i], IMRECT);
264
265 /*
266 * Temporary deletion until above fix tested - GAB 18 Feb 2003
267 * If test works delete permanently.
268 *
269 stack_push((void*) ims[i],IMRECT,im_free);
270 imcalc_draw(imcalc_tv_get());
271 */
272 }
273 pvector_free(ims,seq_start-1);
274 dvector_free(factor, seq_start-1);
275
276 return imlistptr;
277 }
278
279
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.