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

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

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