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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedOctants_csfcount12.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    :
 37  * Date    :
 38  * Version :
 39  * CVS Id  :
 40  *
 41  * Notes :
 42  *
 43  *********
 44 */
 45 
 46 #include "tlmedOctants_csfcount12.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #    include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <float.h>
 54 #include <limits.h>
 55 #include <tina/sys/sysDef.h>
 56 #include <tina/sys/sysPro.h>
 57 #include <tina/math/mathDef.h>
 58 #include <tina/math/mathPro.h>
 59 #include <tina/image/imgPro.h>
 60 #include <tinatool/draw/drawPro.h>
 61 #include <tinatool/tlmedical/tlmed_CoregPro.h>
 62 #include <tinatool/tlmedical/tlmedOctants_alloc.h>
 63 
 64 /* Update of csf_count12 to use box center points rather than box boundaries: store these center points in a list */
 65 /*                                              PAB 10 / 9 / 2002                                                 */
 66 /* Dimensionality variables are as follows: iimptrxx are csf extents in the original data space                   */
 67 /*                                          nimptrxx are csf extents in the normal brain space                    */
 68 /*                                          timptrxx are the coreg tool roi dimensions                            */
 69 /*                                          cimptrxx are the measured normal brain csf extents in the normal      */
 70 /*                                          brain space                                                           */
 71 
 72 static void ***imptrs=NULL;
 73 int nimptrlx, nimptrux, nimptrly, nimptruy, nimptrlz, nimptruz;
 74 
 75 /* Normal brain csf extents in normal brain space for fixed centre points: */
 76 
 77 int cimptrlx = 55;
 78 int cimptrux = 198;
 79 int cimptrly = 57;
 80 int cimptruy = 166;
 81 int cimptrlz = 1;
 82 int cimptruz = 45;
 83 
 84 
 85 static double distance_to_center_point(double *data_point, double *center_point)
 86 {
 87         int i;
 88         double dist=0;
 89 
 90         for(i=0; i<3; i++)
 91         {
 92                 dist += sqr(data_point[i] - center_point[i]);
 93         }
 94         return(dist);
 95 }
 96 
 97 static Match *nearest_center_point(Vec3 post, List *center_point_list, int *side_flag)
 98 {
 99         double *data_point, *center_point;
100         double dist, min_dist = FLT_MAX, x_center;
101         List *plist;
102         Match *nptr, *min_ptr=NULL;
103         int i;
104 
105         if(center_point_list==NULL) return(NULL);
106 
107         data_point = dvector_alloc(0, 3);
108         for(i=0; i<3; i++) data_point[i] = post.el[i];
109 
110         x_center = (double)(cimptrlx+cimptrux)/2.0;
111 
112         if(data_point[0]>(x_center+0.5))
113         {
114                 data_point[0] = (double)(2*x_center - data_point[0]);
115                 *side_flag = 2.0;
116         }       
117         else
118         {
119                 *side_flag = 1.0;
120         }
121         
122         for(plist = center_point_list; plist!=NULL; plist = plist->next)
123         {
124                 nptr = (Match *)plist->to;
125                 center_point = (double *)nptr->to1;
126                 dist = distance_to_center_point(data_point, center_point);
127                 if(dist<=min_dist)
128                 {
129                         min_dist = dist;
130                         min_ptr = nptr;
131                 }
132         }
133 
134         fvector_free(data_point, 0);
135         return(min_ptr);
136 }
137 
138 
139 List *csf_count12_octants_proc(float *brain_dim, List *center_point_list)
140 {
141         Vec3 posi,post;
142         int i,j,k;
143         Imrect *xcoreg_im;
144         double scale;
145         int imptrlx, imptrux,imptrly, imptruy, imptrlz, imptruz;
146         Vec3 iimptrlx, iimptrux, iimptrly, iimptruy, iimptrlz, iimptruz, coreg_c;
147         float timptrlx, timptrux, timptrly, timptruy, timptrlz, timptruz;
148         double xcentre = 128.0;
149         double ycentre = 128.0;
150         double zcentre = 40.0;
151         double xscale = 1.0;
152         double yscale = 1.0;
153         double zscale = 2.0;
154         double nxscale = 1.0;
155         double nyscale = 1.0;
156         double nzscale = 2.0;
157         float nearest_pixel(void ***rasptrs, Vec3 pos);
158         float xdist, ydist, zdist, *weight_vector;
159         Match *nptr;
160         List *plist;
161         int side_flag=0;
162 
163         imptrs = seq_limits(&imptrlz,&imptruz,&imptrly,&imptruy,&imptrlx,&imptrux);
164         /* init_interp(imptrlx,imptrux, imptrly,imptruy, imptrlz,imptruz);      */              /********* Why is this being called here ? */
165         coreg_scales(&xscale,&yscale,&zscale,&nxscale,&nyscale,&nzscale);
166         coreg_centre(&xcentre,&ycentre,&zcentre,&coreg_c);
167 
168         tv_get_roi(xcoreg_tv_get(), &timptrlz, &timptrly, &timptruz, &timptruy);
169 
170         nimptrlz =  INT_MAX;
171         nimptruz = -INT_MAX;
172         nimptrly =  INT_MAX;
173         nimptruy = -INT_MAX;
174         nimptrlx =  INT_MAX;
175         nimptrux = -INT_MAX;
176 
177         /* loop over reslice coordinates locating maximal extent of CSF */
178 
179         for (k=imptrlz;k<imptruz;k++)
180         {
181                 for (i=imptrly;i<imptruy;i++)
182                 {
183                         for (j=imptrlx;j<imptrux;j++)
184                         {
185                                 posi.el[0] = ((float)j+0.5 );
186                                 posi.el[1] = ((float)i+0.5 );
187                                 posi.el[2] = ((float)k+0.5 );
188 
189                                 if(nearest_pixel(imptrs,posi)>0.0)
190                                 {
191                                         post = coreg_proj((double)posi.el[0],(double)posi.el[1],(double)posi.el[2]);
192                                         if (tina_int(post.el[2])<nimptrlz)
193                                         {
194                                                 nimptrlz = tina_int(post.el[2]);
195                                                 iimptrlz = coreg_bproj(0.0,0.0,post.el[2]);
196                                         }
197                                         if (tina_int(post.el[2])>nimptruz)
198                                         {
199                                                 nimptruz = tina_int(post.el[2]);
200                                                 iimptruz = coreg_bproj(0.0,0.0,post.el[2]);
201                                         }
202                                         if (tina_int(post.el[1])<nimptrly)
203                                         {
204                                                 nimptrly = tina_int(post.el[1]);
205                                                 iimptrly = coreg_bproj(0.0,post.el[1],0.0);
206                                         }
207                                         if (tina_int(post.el[1])>nimptruy)
208                                         {
209                                                 nimptruy = tina_int(post.el[1]);
210                                                 iimptruy = coreg_bproj(0.0,post.el[1],0.0);
211                                         }
212                                         if (tina_int(post.el[0])<nimptrlx)
213                                         {
214                                                 nimptrlx = tina_int(post.el[0]);
215                                                 iimptrlx = coreg_bproj(post.el[0],0.0,0.0);
216                                         }
217                                         if (tina_int(post.el[0])>nimptrux)
218                                         {
219                                                 nimptrux = tina_int(post.el[0]);
220                                                 iimptrux = coreg_bproj(post.el[0],0.0,0.0);
221                                         }
222                                 }
223                         }
224                 }
225         }
226         xcoreg_im = get_xcoreg_im();
227         if (xcoreg_im==NULL) nimptruy = (int)((timptruy-ycentre)*yscale)/nyscale + ycentre;
228 
229         /* If center_point_list is NULL in the call to this function, then the call is from the octants tool and the    */
230         /* center point list is initialised with the fixed normal brain extents above.  If it is not NULL, then the     */
231         /* call is from the GA optimisation (with new center points) and the list should be left alone, but the weights */
232         /* (which carry the csf counts) must be reset to zero                                                           */
233 
234         if(center_point_list==NULL)
235         {
236                 /* center_point_list = center_point_init(cimptrlx, cimptrux, cimptrly, cimptruy, cimptrlz, cimptruz); */ /* Use fixed CSF extents */
237                 center_point_list = center_point_init(nimptrlx, nimptrux, nimptrly, nimptruy, nimptrlz, nimptruz); /* Use extracted CSF extents */
238         }
239         else
240         {
241                 for(plist=center_point_list; plist!=NULL; plist=plist->next)
242                 {
243                         nptr = (Match *)plist->to;
244                         weight_vector = (float *)nptr->to2;
245                         for(i=0; i<2; i++)
246                         {
247                                 weight_vector[i] = 0.0;
248                         }
249                 }
250         }
251 
252         /* now divide up the volume and count the contents */
253 
254         for (k=imptrlz;k<imptruz;k++)
255         {
256                 for (i=imptrly;i<imptruy;i++)
257                 {
258                         for (j=imptrlx;j<imptrux;j++)
259                         {
260                                 posi.el[0] = ((float)j+0.5 );
261                                 posi.el[1] = ((float)i+0.5 );
262                                 posi.el[2] = ((float)k+0.5 );
263 
264                                 if (nearest_pixel(imptrs,posi)!=0.0)
265                                 {
266                                         post = coreg_proj((double)posi.el[0],(double)posi.el[1],(double)posi.el[2]);
267                                         side_flag=0;
268                                         nptr = nearest_center_point(post, center_point_list, &side_flag);
269 
270                                         weight_vector = (float *)nptr->to2;
271                                         if(side_flag==1)
272                                         {
273                                                 weight_vector[0] += nearest_pixel(imptrs, posi);
274                                         }
275                                         else if(side_flag==2)
276                                         {
277                                                 weight_vector[1] += nearest_pixel(imptrs, posi);
278                                         }
279                                 }
280                         }
281                 }
282         }
283 
284         /* now multiply numbers of voxels by voxel sizes (of brain loaded in as 'reslice',
285         since this gives the original voxel sizes now that the air file has been swapped)
286         to get volumes in cubic mm */
287 
288         seq_get_iscales(&xscale, &yscale, &zscale);
289         scale = xscale*yscale*zscale;
290 
291         /* In order to be able to scale these volumes into the standard brain space we need
292         to know the total dimensions of the brain being analyzed in mm - measured as the
293         distance between extreme csf voxels distance measured in the standard frame, scaled
294         back into the frame of data being analyzed, to give absolute distance in mm */
295 
296         xdist = sqrt(xscale*xscale*
297                 (iimptrux.el[0]-iimptrlx.el[0])*(iimptrux.el[0]-iimptrlx.el[0])
298                 + yscale*yscale*
299                 (iimptrux.el[1]-iimptrlx.el[1])*(iimptrux.el[1]-iimptrlx.el[1])
300                 + zscale*zscale*
301                 (iimptrux.el[2]-iimptrlx.el[2])*(iimptrux.el[2]-iimptrlx.el[2]));
302         ydist = sqrt(xscale*xscale*
303                 (iimptruy.el[0]-iimptrly.el[0])*(iimptruy.el[0]-iimptrly.el[0])
304                 + yscale*yscale*
305                 (iimptruy.el[1]-iimptrly.el[1])*(iimptruy.el[1]-iimptrly.el[1])
306                 + zscale*zscale*
307                 (iimptruy.el[2]-iimptrly.el[2])*(iimptruy.el[2]-iimptrly.el[2]));
308         zdist = sqrt(xscale*xscale*
309                 (iimptruz.el[0]-iimptrlz.el[0])*(iimptruz.el[0]-iimptrlz.el[0])
310                 + yscale*yscale*
311                 (iimptruz.el[1]-iimptrlz.el[1])*(iimptruz.el[1]-iimptrlz.el[1])
312                 + zscale*zscale*
313                 (iimptruz.el[2]-iimptrlz.el[2])*(iimptruz.el[2]-iimptrlz.el[2]));
314 
315         if(brain_dim!=NULL)
316         {
317                 brain_dim[0] = xdist;
318                 brain_dim[1] = ydist;
319                 brain_dim[2] = zdist;
320         }
321 
322         for(plist = center_point_list; plist!=NULL; plist=plist->next)
323         {
324                 nptr = (Match *)plist->to;
325                 weight_vector = nptr->to2;
326                 weight_vector[0] *=scale;
327                 weight_vector[1] *=scale;
328         }
329 
330         return center_point_list;
331 }
332 

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