Mercurial > hg > orthanc-stone

comparison Framework/Volumes/ImageBuffer3D.cpp @ 684:7719eb852dd5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
new class: VolumeImageGeometry
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 16 May 2019 16:47:46 +0200
parents dbc1d8bfc68a
children 9a474e90e832
comparison
equal deleted inserted replaced
683:dbc1d8bfc68a 684:7719eb852dd5
116 height_(height), 116 height_(height),
117 depth_(depth), 117 depth_(depth),
118 computeRange_(computeRange), 118 computeRange_(computeRange),
119 hasRange_(false) 119 hasRange_(false)
120 { 120 {
121 LinearAlgebra::AssignVector(voxelDimensions_, 1, 1, 1); 121 geometry_.SetSize(width, height, depth);
122 122
123 LOG(INFO) << "Created a 3D image of size " << width << "x" << height 123 LOG(INFO) << "Created a 3D image of size " << width << "x" << height
124 << "x" << depth << " in " << Orthanc::EnumerationToString(format) 124 << "x" << depth << " in " << Orthanc::EnumerationToString(format)
125 << " (" << (GetEstimatedMemorySize() / (1024ll * 1024ll)) << "MB)"; 125 << " (" << (GetEstimatedMemorySize() / (1024ll * 1024ll)) << "MB)";
126
127 UpdateGeometry();
128 } 126 }
129 127
130 128
131 void ImageBuffer3D::Clear() 129 void ImageBuffer3D::Clear()
132 { 130 {
133 memset(image_.GetBuffer(), 0, image_.GetHeight() * image_.GetPitch()); 131 memset(image_.GetBuffer(), 0, image_.GetHeight() * image_.GetPitch());
134 } 132 }
135 133
136 134
137 void ImageBuffer3D::SetAxialGeometry(const CoordinateSystem3D& geometry) 135
138 { 136
139 axialGeometry_ = geometry; 137 ParallelSlices* ImageBuffer3D::GetGeometry(VolumeProjection projection) const
140 UpdateGeometry(); 138 {
141 } 139 const Vector dimensions = geometry_.GetVoxelDimensions(VolumeProjection_Axial);
142 140 const CoordinateSystem3D& axial = geometry_.GetAxialGeometry();
143 141
144 void ImageBuffer3D::SetVoxelDimensions(double x, 142 std::auto_ptr<ParallelSlices> result(new ParallelSlices);
145 double y, 143
146 double z)
147 {
148 if (x <= 0 ||
149 y <= 0 ||
150 z <= 0)
151 {
152 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
153 }
154 else
155 {
156 LinearAlgebra::AssignVector(voxelDimensions_, x, y, z);
157 UpdateGeometry();
158 }
159 }
160
161
162 Vector ImageBuffer3D::GetVoxelDimensions(VolumeProjection projection) const
163 {
164 switch (projection) 144 switch (projection)
165 { 145 {
166 case VolumeProjection_Axial: 146 case VolumeProjection_Axial:
167 return voxelDimensions_; 147 for (unsigned int z = 0; z < depth_; z++)
168 148 {
169 case VolumeProjection_Coronal: 149 Vector origin = axial.GetOrigin();
170 return LinearAlgebra::CreateVector(voxelDimensions_[0], voxelDimensions_[2], voxelDimensions_[1]); 150 origin += static_cast<double>(z) * dimensions[2] * axial.GetNormal();
171 151
172 case VolumeProjection_Sagittal: 152 result->AddSlice(origin,
173 return LinearAlgebra::CreateVector(voxelDimensions_[1], voxelDimensions_[2], voxelDimensions_[0]); 153 axial.GetAxisX(),
174 154 axial.GetAxisY());
175 default: 155 }
176 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 156 break;
177 } 157
178 } 158 case VolumeProjection_Coronal:
179 159 for (unsigned int y = 0; y < height_; y++)
180 160 {
181 void ImageBuffer3D::GetSliceSize(unsigned int& width, 161 Vector origin = axial.GetOrigin();
182 unsigned int& height, 162 origin += static_cast<double>(y) * dimensions[1] * axial.GetAxisY();
183 VolumeProjection projection) 163 origin += static_cast<double>(depth_ - 1) * dimensions[2] * axial.GetNormal();
184 { 164
185 switch (projection) 165 result->AddSlice(origin,
186 { 166 axial.GetAxisX(),
187 case VolumeProjection_Axial: 167 -axial.GetNormal());
188 width = width_; 168 }
189 height = height_; 169 break;
190 break; 170
191 171 case VolumeProjection_Sagittal:
192 case VolumeProjection_Coronal: 172 for (unsigned int x = 0; x < width_; x++)
193 width = width_; 173 {
194 height = depth_; 174 Vector origin = axial.GetOrigin();
195 break; 175 origin += static_cast<double>(x) * dimensions[0] * axial.GetAxisX();
196 176 origin += static_cast<double>(depth_ - 1) * dimensions[2] * axial.GetNormal();
197 case VolumeProjection_Sagittal: 177
198 width = height_; 178 result->AddSlice(origin,
199 height = depth_; 179 axial.GetAxisY(),
200 break; 180 -axial.GetNormal());
201 181 }
202 default: 182 break;
203 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 183
204 } 184 default:
205 } 185 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
206
207
208 ParallelSlices* ImageBuffer3D::GetGeometry(VolumeProjection projection) const
209 {
210 std::auto_ptr<ParallelSlices> result(new ParallelSlices);
211
212 switch (projection)
213 {
214 case VolumeProjection_Axial:
215 for (unsigned int z = 0; z < depth_; z++)
216 {
217 Vector origin = axialGeometry_.GetOrigin();
218 origin += static_cast<double>(z) * voxelDimensions_[2] * axialGeometry_.GetNormal();
219
220 result->AddSlice(origin,
221 axialGeometry_.GetAxisX(),
222 axialGeometry_.GetAxisY());
223 }
224 break;
225
226 case VolumeProjection_Coronal:
227 for (unsigned int y = 0; y < height_; y++)
228 {
229 Vector origin = axialGeometry_.GetOrigin();
230 origin += static_cast<double>(y) * voxelDimensions_[1] * axialGeometry_.GetAxisY();
231 origin += static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal();
232
233 result->AddSlice(origin,
234 axialGeometry_.GetAxisX(),
235 -axialGeometry_.GetNormal());
236 }
237 break;
238
239 case VolumeProjection_Sagittal:
240 for (unsigned int x = 0; x < width_; x++)
241 {
242 Vector origin = axialGeometry_.GetOrigin();
243 origin += static_cast<double>(x) * voxelDimensions_[0] * axialGeometry_.GetAxisX();
244 origin += static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal();
245
246 result->AddSlice(origin,
247 axialGeometry_.GetAxisY(),
248 -axialGeometry_.GetNormal());
249 }
250 break;
251
252 default:
253 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
254 } 186 }
255 187
256 return result.release(); 188 return result.release();
257 } 189 }
258 190
274 206
275 float sliceMin, sliceMax; 207 float sliceMin, sliceMax;
276 208
277 switch (slice.GetFormat()) 209 switch (slice.GetFormat())
278 { 210 {
279 case Orthanc::PixelFormat_Grayscale8: 211 case Orthanc::PixelFormat_Grayscale8:
280 case Orthanc::PixelFormat_Grayscale16: 212 case Orthanc::PixelFormat_Grayscale16:
281 case Orthanc::PixelFormat_Grayscale32: 213 case Orthanc::PixelFormat_Grayscale32:
282 case Orthanc::PixelFormat_SignedGrayscale16: 214 case Orthanc::PixelFormat_SignedGrayscale16:
283 { 215 {
284 int64_t a, b; 216 int64_t a, b;
285 Orthanc::ImageProcessing::GetMinMaxIntegerValue(a, b, slice); 217 Orthanc::ImageProcessing::GetMinMaxIntegerValue(a, b, slice);
286 sliceMin = static_cast<float>(a); 218 sliceMin = static_cast<float>(a);
287 sliceMax = static_cast<float>(b); 219 sliceMax = static_cast<float>(b);
288 break; 220 break;
289 } 221 }
290 222
291 case Orthanc::PixelFormat_Float32: 223 case Orthanc::PixelFormat_Float32:
292 Orthanc::ImageProcessing::GetMinMaxFloatValue(sliceMin, sliceMax, slice); 224 Orthanc::ImageProcessing::GetMinMaxFloatValue(sliceMin, sliceMax, slice);
293 break; 225 break;
294 226
295 default: 227 default:
296 return; 228 return;
297 } 229 }
298 230
299 if (hasRange_) 231 if (hasRange_)
300 { 232 {
301 minValue_ = std::min(minValue_, sliceMin); 233 minValue_ = std::min(minValue_, sliceMin);
365 case VolumeProjection_Sagittal: 297 case VolumeProjection_Sagittal:
366 sagittal_.reset(that.ExtractSagittalSlice(slice)); 298 sagittal_.reset(that.ExtractSagittalSlice(slice));
367 sagittal_->GetReadOnlyAccessor(accessor_); 299 sagittal_->GetReadOnlyAccessor(accessor_);
368 break; 300 break;
369 301
370 default: 302 default:
371 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 303 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
372 } 304 }
373 } 305 }
374 306
375 307
376 void ImageBuffer3D::SliceWriter::Flush() 308 void ImageBuffer3D::SliceWriter::Flush()
409 case VolumeProjection_Sagittal: 341 case VolumeProjection_Sagittal:
410 sagittal_.reset(that.ExtractSagittalSlice(slice)); 342 sagittal_.reset(that.ExtractSagittalSlice(slice));
411 sagittal_->GetWriteableAccessor(accessor_); 343 sagittal_->GetWriteableAccessor(accessor_);
412 break; 344 break;
413 345
414 default: 346 default:
415 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 347 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
416 } 348 }
417 } 349 }
418 350
419 351
420 uint8_t ImageBuffer3D::GetVoxelGrayscale8(unsigned int x, 352 uint8_t ImageBuffer3D::GetVoxelGrayscale8(unsigned int x,
455 } 387 }
456 388
457 const void* p = image_.GetConstRow(y + height_ * (depth_ - 1 - z)); 389 const void* p = image_.GetConstRow(y + height_ * (depth_ - 1 - z));
458 return reinterpret_cast<const uint16_t*>(p) [x]; 390 return reinterpret_cast<const uint16_t*>(p) [x];
459 } 391 }
460
461
462 void ImageBuffer3D::UpdateGeometry()
463 {
464 Vector p = (axialGeometry_.GetOrigin() +
465 static_cast<double>(depth_ - 1) * voxelDimensions_[2] * axialGeometry_.GetNormal());
466
467 coronalGeometry_ = CoordinateSystem3D(p,
468 axialGeometry_.GetAxisX(),
469 -axialGeometry_.GetNormal());
470
471 sagittalGeometry_ = CoordinateSystem3D(p,
472 axialGeometry_.GetAxisY(),
473 -axialGeometry_.GetNormal());
474
475 Vector origin = (
476 axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],
477 -0.5 * voxelDimensions_[1]) -
478 0.5 * voxelDimensions_[2] * axialGeometry_.GetNormal());
479
480 Vector scaling = (
481 axialGeometry_.GetAxisX() * voxelDimensions_[0] * static_cast<double>(width_) +
482 axialGeometry_.GetAxisY() * voxelDimensions_[1] * static_cast<double>(height_) +
483 axialGeometry_.GetNormal() * voxelDimensions_[2] * static_cast<double>(depth_));
484
485 transform_ = LinearAlgebra::Product(
486 GeometryToolbox::CreateTranslationMatrix(origin[0], origin[1], origin[2]),
487 GeometryToolbox::CreateScalingMatrix(scaling[0], scaling[1], scaling[2]));
488
489 LinearAlgebra::InvertMatrix(transformInverse_, transform_);
490 }
491
492
493 Vector ImageBuffer3D::GetCoordinates(float x,
494 float y,
495 float z) const
496 {
497 Vector p = LinearAlgebra::Product(transform_, LinearAlgebra::CreateVector(x, y, z, 1));
498
499 assert(LinearAlgebra::IsNear(p[3], 1)); // Affine transform, no perspective effect
500
501 // Back to non-homogeneous coordinates
502 return LinearAlgebra::CreateVector(p[0], p[1], p[2]);
503 }
504
505
506 bool ImageBuffer3D::DetectProjection(VolumeProjection& projection,
507 const CoordinateSystem3D& plane) const
508 {
509 if (GeometryToolbox::IsParallel(plane.GetNormal(),
510 axialGeometry_.GetNormal()))
511 {
512 projection = VolumeProjection_Axial;
513 return true;
514 }
515
516 {
517 if (GeometryToolbox::IsParallel(plane.GetNormal(),
518 coronalGeometry_.GetNormal()))
519 {
520 projection = VolumeProjection_Coronal;
521 return true;
522 }
523 }
524
525 {
526 if (GeometryToolbox::IsParallel(plane.GetNormal(),
527 sagittalGeometry_.GetNormal()))
528 {
529 projection = VolumeProjection_Sagittal;
530 return true;
531 }
532 }
533
534 return false;
535 }
536 } 392 }