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

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

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