00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 namespace PolyVox
00025 {
00026 template< template<typename> class VolumeType, typename VoxelType>
00027 Vector3DFloat computeCentralDifferenceGradient(const typename VolumeType<VoxelType>::Sampler& volIter)
00028 {
00029
00030 VoxelType voxel1nx = volIter.peekVoxel1nx0py0pz() > 0 ? 1: 0;
00031 VoxelType voxel1px = volIter.peekVoxel1px0py0pz() > 0 ? 1: 0;
00032
00033 VoxelType voxel1ny = volIter.peekVoxel0px1ny0pz() > 0 ? 1: 0;
00034 VoxelType voxel1py = volIter.peekVoxel0px1py0pz() > 0 ? 1: 0;
00035
00036 VoxelType voxel1nz = volIter.peekVoxel0px0py1nz() > 0 ? 1: 0;
00037 VoxelType voxel1pz = volIter.peekVoxel0px0py1pz() > 0 ? 1: 0;
00038
00039 return Vector3DFloat
00040 (
00041 static_cast<float>(voxel1nx) - static_cast<float>(voxel1px),
00042 static_cast<float>(voxel1ny) - static_cast<float>(voxel1py),
00043 static_cast<float>(voxel1nz) - static_cast<float>(voxel1pz)
00044 );
00045 }
00046
00047 template< template<typename> class VolumeType, typename VoxelType>
00048 Vector3DFloat computeDecimatedCentralDifferenceGradient(const typename VolumeType<VoxelType>::Sampler& volIter)
00049 {
00050 const int32_t x = volIter.getPosX();
00051 const int32_t y = volIter.getPosY();
00052 const int32_t z = volIter.getPosZ();
00053
00054
00055 VoxelType voxel1nx = volIter.getVoxelAt(x-2, y ,z ) > 0 ? 1: 0;
00056 VoxelType voxel1px = volIter.getVoxelAt(x-2, y ,z ) > 0 ? 1: 0;
00057
00058 VoxelType voxel1ny = volIter.getVoxelAt(x , y-2,z ) > 0 ? 1: 0;
00059 VoxelType voxel1py = volIter.getVoxelAt(x , y-2,z ) > 0 ? 1: 0;
00060
00061 VoxelType voxel1nz = volIter.getVoxelAt(x , y ,z-2) > 0 ? 1: 0;
00062 VoxelType voxel1pz = volIter.getVoxelAt(x , y ,z-2) > 0 ? 1: 0;
00063
00064 return Vector3DFloat
00065 (
00066 static_cast<float>(voxel1nx) - static_cast<float>(voxel1px),
00067 static_cast<float>(voxel1ny) - static_cast<float>(voxel1py),
00068 static_cast<float>(voxel1nz) - static_cast<float>(voxel1pz)
00069 );
00070 }
00071
00072 template< template<typename> class VolumeType, typename VoxelType>
00073 Vector3DFloat computeSmoothCentralDifferenceGradient(typename VolumeType<VoxelType>::Sampler& volIter)
00074 {
00075 int32_t initialX = volIter.getPosX();
00076 int32_t initialY = volIter.getPosY();
00077 int32_t initialZ = volIter.getPosZ();
00078
00079
00080 volIter.setPosition(initialX-1, initialY, initialZ);
00081 float voxel1nx = computeSmoothedVoxel(volIter);
00082 volIter.setPosition(initialX+1, initialY, initialZ);
00083 float voxel1px = computeSmoothedVoxel(volIter);
00084
00085 volIter.setPosition(initialX, initialY-1, initialZ);
00086 float voxel1ny = computeSmoothedVoxel(volIter);
00087 volIter.setPosition(initialX, initialY+1, initialZ);
00088 float voxel1py = computeSmoothedVoxel(volIter);
00089
00090 volIter.setPosition(initialX, initialY, initialZ-1);
00091 float voxel1nz = computeSmoothedVoxel(volIter);
00092 volIter.setPosition(initialX, initialY, initialZ+1);
00093 float voxel1pz = computeSmoothedVoxel(volIter);
00094
00095 return Vector3DFloat
00096 (
00097 voxel1nx - voxel1px,
00098 voxel1ny - voxel1py,
00099 voxel1nz - voxel1pz
00100 );
00101 }
00102
00103 template< template<typename> class VolumeType, typename VoxelType>
00104 Vector3DFloat computeSobelGradient(const typename VolumeType<VoxelType>::Sampler& volIter)
00105 {
00106 static const int weights[3][3][3] = { { {2,3,2}, {3,6,3}, {2,3,2} }, {
00107 {3,6,3}, {6,0,6}, {3,6,3} }, { {2,3,2}, {3,6,3}, {2,3,2} } };
00108
00109 const VoxelType pVoxel1nx1ny1nz = volIter.peekVoxel1nx1ny1nz() > 0 ? 1: 0;
00110 const VoxelType pVoxel1nx1ny0pz = volIter.peekVoxel1nx1ny0pz() > 0 ? 1: 0;
00111 const VoxelType pVoxel1nx1ny1pz = volIter.peekVoxel1nx1ny1pz() > 0 ? 1: 0;
00112 const VoxelType pVoxel1nx0py1nz = volIter.peekVoxel1nx0py1nz() > 0 ? 1: 0;
00113 const VoxelType pVoxel1nx0py0pz = volIter.peekVoxel1nx0py0pz() > 0 ? 1: 0;
00114 const VoxelType pVoxel1nx0py1pz = volIter.peekVoxel1nx0py1pz() > 0 ? 1: 0;
00115 const VoxelType pVoxel1nx1py1nz = volIter.peekVoxel1nx1py1nz() > 0 ? 1: 0;
00116 const VoxelType pVoxel1nx1py0pz = volIter.peekVoxel1nx1py0pz() > 0 ? 1: 0;
00117 const VoxelType pVoxel1nx1py1pz = volIter.peekVoxel1nx1py1pz() > 0 ? 1: 0;
00118
00119 const VoxelType pVoxel0px1ny1nz = volIter.peekVoxel0px1ny1nz() > 0 ? 1: 0;
00120 const VoxelType pVoxel0px1ny0pz = volIter.peekVoxel0px1ny0pz() > 0 ? 1: 0;
00121 const VoxelType pVoxel0px1ny1pz = volIter.peekVoxel0px1ny1pz() > 0 ? 1: 0;
00122 const VoxelType pVoxel0px0py1nz = volIter.peekVoxel0px0py1nz() > 0 ? 1: 0;
00123
00124 const VoxelType pVoxel0px0py1pz = volIter.peekVoxel0px0py1pz() > 0 ? 1: 0;
00125 const VoxelType pVoxel0px1py1nz = volIter.peekVoxel0px1py1nz() > 0 ? 1: 0;
00126 const VoxelType pVoxel0px1py0pz = volIter.peekVoxel0px1py0pz() > 0 ? 1: 0;
00127 const VoxelType pVoxel0px1py1pz = volIter.peekVoxel0px1py1pz() > 0 ? 1: 0;
00128
00129 const VoxelType pVoxel1px1ny1nz = volIter.peekVoxel1px1ny1nz() > 0 ? 1: 0;
00130 const VoxelType pVoxel1px1ny0pz = volIter.peekVoxel1px1ny0pz() > 0 ? 1: 0;
00131 const VoxelType pVoxel1px1ny1pz = volIter.peekVoxel1px1ny1pz() > 0 ? 1: 0;
00132 const VoxelType pVoxel1px0py1nz = volIter.peekVoxel1px0py1nz() > 0 ? 1: 0;
00133 const VoxelType pVoxel1px0py0pz = volIter.peekVoxel1px0py0pz() > 0 ? 1: 0;
00134 const VoxelType pVoxel1px0py1pz = volIter.peekVoxel1px0py1pz() > 0 ? 1: 0;
00135 const VoxelType pVoxel1px1py1nz = volIter.peekVoxel1px1py1nz() > 0 ? 1: 0;
00136 const VoxelType pVoxel1px1py0pz = volIter.peekVoxel1px1py0pz() > 0 ? 1: 0;
00137 const VoxelType pVoxel1px1py1pz = volIter.peekVoxel1px1py1pz() > 0 ? 1: 0;
00138
00139 const int xGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
00140 weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
00141 pVoxel1nx1ny1pz - weights[0][1][0] * pVoxel1nx0py1nz -
00142 weights[1][1][0] * pVoxel1nx0py0pz - weights[2][1][0] *
00143 pVoxel1nx0py1pz - weights[0][2][0] * pVoxel1nx1py1nz -
00144 weights[1][2][0] * pVoxel1nx1py0pz - weights[2][2][0] *
00145 pVoxel1nx1py1pz + weights[0][0][2] * pVoxel1px1ny1nz +
00146 weights[1][0][2] * pVoxel1px1ny0pz + weights[2][0][2] *
00147 pVoxel1px1ny1pz + weights[0][1][2] * pVoxel1px0py1nz +
00148 weights[1][1][2] * pVoxel1px0py0pz + weights[2][1][2] *
00149 pVoxel1px0py1pz + weights[0][2][2] * pVoxel1px1py1nz +
00150 weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
00151 pVoxel1px1py1pz);
00152
00153 const int yGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
00154 weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
00155 pVoxel1nx1ny1pz + weights[0][2][0] * pVoxel1nx1py1nz +
00156 weights[1][2][0] * pVoxel1nx1py0pz + weights[2][2][0] *
00157 pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz -
00158 weights[1][0][1] * pVoxel0px1ny0pz - weights[2][0][1] *
00159 pVoxel0px1ny1pz + weights[0][2][1] * pVoxel0px1py1nz +
00160 weights[1][2][1] * pVoxel0px1py0pz + weights[2][2][1] *
00161 pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz -
00162 weights[1][0][2] * pVoxel1px1ny0pz - weights[2][0][2] *
00163 pVoxel1px1ny1pz + weights[0][2][2] * pVoxel1px1py1nz +
00164 weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
00165 pVoxel1px1py1pz);
00166
00167 const int zGrad(- weights[0][0][0] * pVoxel1nx1ny1nz +
00168 weights[2][0][0] * pVoxel1nx1ny1pz - weights[0][1][0] *
00169 pVoxel1nx0py1nz + weights[2][1][0] * pVoxel1nx0py1pz -
00170 weights[0][2][0] * pVoxel1nx1py1nz + weights[2][2][0] *
00171 pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz +
00172 weights[2][0][1] * pVoxel0px1ny1pz - weights[0][1][1] *
00173 pVoxel0px0py1nz + weights[2][1][1] * pVoxel0px0py1pz -
00174 weights[0][2][1] * pVoxel0px1py1nz + weights[2][2][1] *
00175 pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz +
00176 weights[2][0][2] * pVoxel1px1ny1pz - weights[0][1][2] *
00177 pVoxel1px0py1nz + weights[2][1][2] * pVoxel1px0py1pz -
00178 weights[0][2][2] * pVoxel1px1py1nz + weights[2][2][2] *
00179 pVoxel1px1py1pz);
00180
00181
00182
00183 return Vector3DFloat(static_cast<float>(-xGrad),static_cast<float>(-yGrad),static_cast<float>(-zGrad));
00184 }
00185
00186 template< template<typename> class VolumeType, typename VoxelType>
00187 Vector3DFloat computeSmoothSobelGradient(typename VolumeType<VoxelType>::Sampler& volIter)
00188 {
00189 static const int weights[3][3][3] = { { {2,3,2}, {3,6,3}, {2,3,2} }, {
00190 {3,6,3}, {6,0,6}, {3,6,3} }, { {2,3,2}, {3,6,3}, {2,3,2} } };
00191
00192 int32_t initialX = volIter.getPosX();
00193 int32_t initialY = volIter.getPosY();
00194 int32_t initialZ = volIter.getPosZ();
00195
00196 volIter.setPosition(initialX-1, initialY-1, initialZ-1); const float pVoxel1nx1ny1nz = computeSmoothedVoxel(volIter);
00197 volIter.setPosition(initialX-1, initialY-1, initialZ ); const float pVoxel1nx1ny0pz = computeSmoothedVoxel(volIter);
00198 volIter.setPosition(initialX-1, initialY-1, initialZ+1); const float pVoxel1nx1ny1pz = computeSmoothedVoxel(volIter);
00199 volIter.setPosition(initialX-1, initialY , initialZ-1); const float pVoxel1nx0py1nz = computeSmoothedVoxel(volIter);
00200 volIter.setPosition(initialX-1, initialY , initialZ ); const float pVoxel1nx0py0pz = computeSmoothedVoxel(volIter);
00201 volIter.setPosition(initialX-1, initialY , initialZ+1); const float pVoxel1nx0py1pz = computeSmoothedVoxel(volIter);
00202 volIter.setPosition(initialX-1, initialY+1, initialZ-1); const float pVoxel1nx1py1nz = computeSmoothedVoxel(volIter);
00203 volIter.setPosition(initialX-1, initialY+1, initialZ ); const float pVoxel1nx1py0pz = computeSmoothedVoxel(volIter);
00204 volIter.setPosition(initialX-1, initialY+1, initialZ+1); const float pVoxel1nx1py1pz = computeSmoothedVoxel(volIter);
00205
00206 volIter.setPosition(initialX , initialY-1, initialZ-1); const float pVoxel0px1ny1nz = computeSmoothedVoxel(volIter);
00207 volIter.setPosition(initialX , initialY-1, initialZ ); const float pVoxel0px1ny0pz = computeSmoothedVoxel(volIter);
00208 volIter.setPosition(initialX , initialY-1, initialZ+1); const float pVoxel0px1ny1pz = computeSmoothedVoxel(volIter);
00209 volIter.setPosition(initialX , initialY , initialZ-1); const float pVoxel0px0py1nz = computeSmoothedVoxel(volIter);
00210
00211 volIter.setPosition(initialX , initialY , initialZ+1); const float pVoxel0px0py1pz = computeSmoothedVoxel(volIter);
00212 volIter.setPosition(initialX , initialY+1, initialZ-1); const float pVoxel0px1py1nz = computeSmoothedVoxel(volIter);
00213 volIter.setPosition(initialX , initialY+1, initialZ ); const float pVoxel0px1py0pz = computeSmoothedVoxel(volIter);
00214 volIter.setPosition(initialX , initialY+1, initialZ+1); const float pVoxel0px1py1pz = computeSmoothedVoxel(volIter);
00215
00216 volIter.setPosition(initialX+1, initialY-1, initialZ-1); const float pVoxel1px1ny1nz = computeSmoothedVoxel(volIter);
00217 volIter.setPosition(initialX+1, initialY-1, initialZ ); const float pVoxel1px1ny0pz = computeSmoothedVoxel(volIter);
00218 volIter.setPosition(initialX+1, initialY-1, initialZ+1); const float pVoxel1px1ny1pz = computeSmoothedVoxel(volIter);
00219 volIter.setPosition(initialX+1, initialY , initialZ-1); const float pVoxel1px0py1nz = computeSmoothedVoxel(volIter);
00220 volIter.setPosition(initialX+1, initialY , initialZ ); const float pVoxel1px0py0pz = computeSmoothedVoxel(volIter);
00221 volIter.setPosition(initialX+1, initialY , initialZ+1); const float pVoxel1px0py1pz = computeSmoothedVoxel(volIter);
00222 volIter.setPosition(initialX+1, initialY+1, initialZ-1); const float pVoxel1px1py1nz = computeSmoothedVoxel(volIter);
00223 volIter.setPosition(initialX+1, initialY+1, initialZ ); const float pVoxel1px1py0pz = computeSmoothedVoxel(volIter);
00224 volIter.setPosition(initialX+1, initialY+1, initialZ+1); const float pVoxel1px1py1pz = computeSmoothedVoxel(volIter);
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 const float xGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
00257 weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
00258 pVoxel1nx1ny1pz - weights[0][1][0] * pVoxel1nx0py1nz -
00259 weights[1][1][0] * pVoxel1nx0py0pz - weights[2][1][0] *
00260 pVoxel1nx0py1pz - weights[0][2][0] * pVoxel1nx1py1nz -
00261 weights[1][2][0] * pVoxel1nx1py0pz - weights[2][2][0] *
00262 pVoxel1nx1py1pz + weights[0][0][2] * pVoxel1px1ny1nz +
00263 weights[1][0][2] * pVoxel1px1ny0pz + weights[2][0][2] *
00264 pVoxel1px1ny1pz + weights[0][1][2] * pVoxel1px0py1nz +
00265 weights[1][1][2] * pVoxel1px0py0pz + weights[2][1][2] *
00266 pVoxel1px0py1pz + weights[0][2][2] * pVoxel1px1py1nz +
00267 weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
00268 pVoxel1px1py1pz);
00269
00270 const float yGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
00271 weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
00272 pVoxel1nx1ny1pz + weights[0][2][0] * pVoxel1nx1py1nz +
00273 weights[1][2][0] * pVoxel1nx1py0pz + weights[2][2][0] *
00274 pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz -
00275 weights[1][0][1] * pVoxel0px1ny0pz - weights[2][0][1] *
00276 pVoxel0px1ny1pz + weights[0][2][1] * pVoxel0px1py1nz +
00277 weights[1][2][1] * pVoxel0px1py0pz + weights[2][2][1] *
00278 pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz -
00279 weights[1][0][2] * pVoxel1px1ny0pz - weights[2][0][2] *
00280 pVoxel1px1ny1pz + weights[0][2][2] * pVoxel1px1py1nz +
00281 weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
00282 pVoxel1px1py1pz);
00283
00284 const float zGrad(- weights[0][0][0] * pVoxel1nx1ny1nz +
00285 weights[2][0][0] * pVoxel1nx1ny1pz - weights[0][1][0] *
00286 pVoxel1nx0py1nz + weights[2][1][0] * pVoxel1nx0py1pz -
00287 weights[0][2][0] * pVoxel1nx1py1nz + weights[2][2][0] *
00288 pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz +
00289 weights[2][0][1] * pVoxel0px1ny1pz - weights[0][1][1] *
00290 pVoxel0px0py1nz + weights[2][1][1] * pVoxel0px0py1pz -
00291 weights[0][2][1] * pVoxel0px1py1nz + weights[2][2][1] *
00292 pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz +
00293 weights[2][0][2] * pVoxel1px1ny1pz - weights[0][1][2] *
00294 pVoxel1px0py1nz + weights[2][1][2] * pVoxel1px0py1pz -
00295 weights[0][2][2] * pVoxel1px1py1nz + weights[2][2][2] *
00296 pVoxel1px1py1pz);
00297
00298
00299
00300 return Vector3DFloat(-xGrad,-yGrad,-zGrad);
00301 }
00302 }