aimsalgo  5.0.5
Neuroimaging image processing
curv3Diso.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef AIMS_MATH_CURV3DISO_H
36 #define AIMS_MATH_CURV3DISO_H
37 
38 #include <aims/math/mathelem.h>
39 #include <aims/data/fastAllocationData.h>
40 
43 enum AimsCurvatureType
44 {
45  AIMS_GAUSSIAN_CURVATURE,
46  AIMS_MEAN_CURVATURE
47 };
48 
49 
52 template <class T> inline
53 AimsData<float>
54 AimsIsoIntensityCurvature3D(const AimsData<T> &vol,
55  AimsCurvatureType type = AIMS_MEAN_CURVATURE);
56 
57 
58 template <class T>
59 AimsData<float>
60 AimsIsoIntensityCurvature3D(const AimsData<T> &vol,AimsCurvatureType type)
61 { ASSERT(vol.dimT()==1);
62  aims::AimsFastAllocationData<float> curv(vol.dimX(),vol.dimY(),vol.dimZ(),1,vol.borderWidth());
63  curv.fillBorder(0);
64  curv = 0.0;
65  curv.setSizeX(vol.sizeX());
66  curv.setSizeY(vol.sizeY());
67  curv.setSizeZ(vol.sizeZ());
68  float fx,fy,fz,fxx,fyy,fzz,fxy,fxz,fyz;
69 
70  switch (type)
71  {
72  case AIMS_GAUSSIAN_CURVATURE:
73  for (int z=2;z<vol.dimZ()-2;z++)
74  for (int y=2;y<vol.dimY()-2;y++)
75  for (int x=2;x<vol.dimX()-2;x++)
76  { if (vol(x-1,y,z)!=vol(x+1,y,z) ||
77  vol(x,y-1,z)!=vol(x,y+1,z) ||
78  vol(x,y,z-1)!=vol(x,y,z+1) )
79  { fx = ((float)vol(x+1,y,z) - (float)vol(x-1,y,z) ) / 2;
80  fy = ((float)vol(x,y+1,z) - (float)vol(x,y-1,z) ) / 2;
81  fz = ((float)vol(x,y,z+1) - (float)vol(x,y,z-1) ) / 2;
82  fxx = ((float)vol(x+2,y,z) - 2 * (float)vol(x,y,z) +
83  (float)vol(x-2,y,z) ) / 4;
84  fyy = ((float)vol(x,y+2,z) - 2 * (float)vol(x,y,z) +
85  (float)vol(x,y-2,z) ) / 4;
86  fzz = ((float)vol(x,y,z+2) - 2 * (float)vol(x,y,z) +
87  (float)vol(x,y,z-2) ) / 4;
88  fxy = ((float)vol(x+1,y+1,z) -
89  (float)vol(x-1,y+1,z) -
90  (float)vol(x+1,y-1,z) +
91  (float)vol(x-1,y-1,z) ) / 4;
92  fxz = ((float)vol(x+1,y,z+1) -
93  (float)vol(x-1,y,z+1) -
94  (float)vol(x+1,y,z-1) +
95  (float)vol(x-1,y,z-1) ) / 4;
96  fyz = ((float)vol(x,y+1,z+1) -
97  (float)vol(x,y-1,z+1) -
98  (float)vol(x,y+1,z-1) +
99  (float)vol(x,y-1,z-1) ) / 4;
100  curv(x,y,z) = ( fx*fx*(fyy*fzz-fyz*fyz) +
101  2*fy*fz*(fxz*fxy-fxx*fyz) +
102  fy*fy*(fxx*fzz-fxz*fxz) +
103  2*fx*fz*(fyz*fxy-fyy*fxz) +
104  fz*fz*(fxx*fyy-fxy*fxy) +
105  2*fx*fy*(fxz*fyz-fzz*fxy) ) /
106  ((fx*fx+fy*fy+fz*fz)*(fx*fx+fy*fy+fz*fz));
107  }
108  else curv(x,y,z) = 0;
109  }
110  break;
111  case AIMS_MEAN_CURVATURE:
112  for (int z=2;z<vol.dimZ()-2;z++)
113  for (int y=2;y<vol.dimY()-2;y++)
114  for (int x=2;x<vol.dimX()-2;x++)
115  { if (vol(x-1,y,z)!=vol(x+1,y,z) ||
116  vol(x,y-1,z)!=vol(x,y+1,z) ||
117  vol(x,y,z-1)!=vol(x,y,z+1) )
118  { fx = ((float)vol(x+1,y,z) - (float)vol(x-1,y,z) ) / 2;
119  fy = ((float)vol(x,y+1,z) - (float)vol(x,y-1,z) ) / 2;
120  fz = ((float)vol(x,y,z+1) - (float)vol(x,y,z-1) ) / 2;
121  fxx = ((float)vol(x+2,y,z) - 2 * (float)vol(x,y,z) +
122  (float)vol(x-2,y,z) ) / 4;
123  fyy = ((float)vol(x,y+2,z) - 2 * (float)vol(x,y,z) +
124  (float)vol(x,y-2,z) ) / 4;
125  fzz = ((float)vol(x,y,z+2) - 2 * (float)vol(x,y,z) +
126  (float)vol(x,y,z-2) ) / 4;
127  fxy = ((float)vol(x+1,y+1,z) -
128  (float)vol(x-1,y+1,z) -
129  (float)vol(x+1,y-1,z) +
130  (float)vol(x-1,y-1,z) ) / 4;
131  fxz = ((float)vol(x+1,y,z+1) -
132  (float)vol(x-1,y,z+1) -
133  (float)vol(x+1,y,z-1) +
134  (float)vol(x-1,y,z-1) ) / 4;
135  fyz = ((float)vol(x,y+1,z+1) -
136  (float)vol(x,y-1,z+1) -
137  (float)vol(x,y+1,z-1) +
138  (float)vol(x,y-1,z-1) ) / 4;
139  curv(x,y,z) = ( fx*fx*(fyy+fzz) - 2*fy*fz*fyz +
140  fy*fy*(fxx+fzz) - 2*fx*fz*fxz +
141  fz*fz*(fxx+fyy) - 2*fx*fy*fxy ) /
142  ( 2 * cube( sqrt(fx*fx + fy*fy + fz*fz) ) );
143  }
144  else curv(x,y,z) = 0;
145  }
146  break;
147  }
148  return curv;
149 }
150 
151 
152 #endif
AimsIsoIntensityCurvature3D
AimsData< float > AimsIsoIntensityCurvature3D(const AimsData< T > &vol, AimsCurvatureType type=AIMS_MEAN_CURVATURE)
3D curvature functions on an intensity image f(x,y,z) = I
Definition: curv3Diso.h:60
AimsData::dimZ
int dimZ() const
AimsData::sizeZ
float sizeZ() const
AimsData::dimY
int dimY() const
AimsData::sizeX
float sizeX() const
cube
T cube(const T &val)
AimsCurvatureType
AimsCurvatureType
The different 3D curvature types.
Definition: curv3Diso.h:43
AimsData::fillBorder
void fillBorder(const T &val)
AimsData::sizeY
float sizeY() const
AimsData::dimT
int dimT() const
ASSERT
#define ASSERT(EX)
AimsData::dimX
int dimX() const
aims::Border::borderWidth
int borderWidth() const