From 5161efae5c49a20be46a49da1ad35d41b0318ad8 Mon Sep 17 00:00:00 2001
From: Andreas Schenk <andreas_schenk@hms.harvard.edu>
Date: Sun, 11 Nov 2012 16:08:51 -0500
Subject: [PATCH] refactored ValueHolder, ImageStateImpl and FFT
Changed in detail:
------------------
- removed DEPRECATED and unused code (slice and row pointers) in value holder
- switched from new allocatd array to std::vector
- removed explicit data initialization loop
- using size_t instead of unsigned int for memory allocation
- removed padding code
- removed try catch(...) statements
- added overflow check for big volumes
- optimized 3D index to linear index calculation
- removed physical extent for spatial domain (identical to logical extend not that padding is gone)
- switched FFT to out of place transforms
Benefits
--------
- memory for large images is correctly allocated
- if the memory can't be allocated an exception is thrown instead of crashing the program
- speed increase for FFT ca. 12% and for algorithms using direct value access (e.g. Stat) ca. 7%
Issues
-----
- FFTs are not out of place, leading to increased memory usage for c2c DFT backtransform and
hc2r FFT/DFT backtransform (3x image size instead of 2x image size). If necessary it would
be possible to add a fallback for these cases (i.e. catch bad_alloc for the tmp state) and
do them in place in low memory conditions. Padding for the in-place transform could be added
to the std::vector either at construction time or otherwise using reserve().
---
modules/img/alg/src/fft.cc | 162 ++++-------
modules/img/alg/src/fractional_shift.hh | 2 +-
modules/img/alg/src/histogram.cc | 2 +-
modules/img/alg/src/transcendentals.hh | 12 +-
modules/img/alg/tests/test_fft.cc | 2 +-
.../base/src/image_state/image_state_base.hh | 2 +-
.../base/src/image_state/image_state_impl.cc | 12 +-
.../base/src/image_state/image_state_impl.hh | 19 +-
.../image_state/image_state_spatial_domain.hh | 27 +-
.../img/base/src/image_state/value_holder.cc | 265 +++++-------------
.../img/base/src/image_state/value_holder.hh | 132 ++-------
modules/img/base/tests/test_domains.cc | 25 +-
modules/img/base/tests/test_value_holder.cc | 17 +-
modules/io/src/img/map_io_mrc_handler.cc | 2 +-
14 files changed, 181 insertions(+), 500 deletions(-)
diff --git a/modules/img/alg/src/fft.cc b/modules/img/alg/src/fft.cc
index ddc137d1a..955f21834 100644
--- a/modules/img/alg/src/fft.cc
+++ b/modules/img/alg/src/fft.cc
@@ -78,39 +78,12 @@ template <>
ImageStateBasePtr FFTFnc::VisitState<Real,SpatialDomain>(const RealSpatialImageState& in_state) const
{
int rank,n[3]; // for actual fftw call
- size_t block_count,src_size,dst_size; // padding parameters
Size in_size = in_state.GetExtent().GetSize();
-
- switch (in_size.GetDim()) {
- case 1:
- rank = 1;
- n[0] = in_size.GetWidth();
- block_count=1;
- src_size=n[0]; // real values
- dst_size=half_plus_one(n[0]); // complex values
- break;
- case 2:
- rank = 2;
- n[0] = in_size.GetWidth();
- n[1] = in_size.GetHeight();
- block_count=n[0];
- src_size=n[1]; // real values
- dst_size=half_plus_one(n[1]); // complex values
- break;
- case 3:
- rank = 3;
- n[0] = in_size.GetWidth();
- n[1] = in_size.GetHeight();
- n[2] = in_size.GetDepth();
- block_count=n[0]*n[1];
- src_size=n[2]; // real values
- dst_size=half_plus_one(n[2]); // complex values
- break;
- default:
- throw(FFTException("unexpected dimension in FFT R2C"));
- }
-
+ n[0] = in_size.GetWidth();
+ n[1] = in_size.GetHeight();
+ n[2] = in_size.GetDepth();
+ rank=in_size.GetDim();
PixelSampling ps=in_state.GetSampling();
ps.SetDomain(FREQUENCY);
@@ -121,25 +94,10 @@ ImageStateBasePtr FFTFnc::VisitState<Real,SpatialDomain>(const RealSpatialImageS
out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
- OST_FFTW_fftw_complex* fftw_out =
-reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
-
- /*
- Since the input array is not preserved in fftw3, a copy is made, and
- then an in-place transform is performed. The real array must be padded
- correctly for this to work.
- */
- Real* fftw_in = reinterpret_cast<Real*>(fftw_out);
- Real* in_ptr = reinterpret_cast<Real*>(in_state.Data().GetData());
-
- for(size_t i=0;i<block_count;i++) {
- std::copy(&in_ptr[i*src_size],&in_ptr[(i+1)*src_size],&fftw_in[i*2*dst_size]);
- }
+ OST_FFTW_fftw_complex* fftw_out = reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
+ Real* fftw_in = const_cast<Real*>(in_state.Data().GetData());
OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
- OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft_r2c(rank,n,
- fftw_in,fftw_out,
- FFTW_ESTIMATE);
-
+ OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft_r2c(rank,n,fftw_in,fftw_out,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
OST_FFTW_fftw_execute(plan);
OST_FFTW_fftw_destroy_plan(plan);
@@ -153,9 +111,11 @@ template <>
ImageStateBasePtr FFTFnc::VisitState<Complex,HalfFrequencyDomain>(const ComplexHalfFrequencyImageState& in_state) const
{
int rank,n[3]; // for actual fftw call
- size_t block_count,src_size,dst_size; // padding parameters
Size in_size=in_state.GetLogicalExtent().GetSize();
+ n[0]=in_size.GetWidth();
+ n[1]=in_size.GetHeight();
+ n[2]=in_size.GetDepth();
// copy
ComplexHalfFrequencyImageState tmp_state(in_state);
// correct phase origin
@@ -164,31 +124,16 @@ ImageStateBasePtr FFTFnc::VisitState<Complex,HalfFrequencyDomain>(const ComplexH
switch (in_size.GetDim()) {
case 1:
rank=1;
- n[0]=in_size.GetWidth();
- block_count=0; // no un-padding for 1D necessary
- src_size=half_plus_one(n[0]); // complex values
- dst_size=n[0]; // real values
break;
case 2:
rank=2;
- n[0]=in_size.GetWidth();
- n[1]=in_size.GetHeight();
- block_count=n[0];
- src_size=half_plus_one(n[1]); // complex values
- dst_size=n[1]; // real values
// complete zero line
for(int i=1;i<half_plus_one(n[0]);++i){
tmp_state.Data().Value(Index(n[0]-i,0,0))=conj(tmp_state.Data().Value(Index(i,0,0)));
}
break;
case 3:
- rank=3;
- n[0]=in_size.GetWidth();
- n[1]=in_size.GetHeight();
- n[2]=in_size.GetDepth();
- block_count=n[0]*n[1];
- src_size=half_plus_one(n[2]); // complex values
- dst_size=n[2]; // real values
+ rank=3;
// complete zero line
for(int i=1;i<half_plus_one(n[0]);++i){
tmp_state.Data().Value(Index(n[0]-i,0,0))=conj(tmp_state.Data().Value(Index(i,0,0)));
@@ -208,9 +153,8 @@ ImageStateBasePtr FFTFnc::VisitState<Complex,HalfFrequencyDomain>(const ComplexH
PixelSampling ps=in_state.GetSampling();
ps.SetDomain(SPATIAL);
- // convert (!) to real spatial image state
Size out_size = in_state.GetLogicalExtent().GetSize();
- boost::shared_ptr<RealSpatialImageState> out_state(new RealSpatialImageState(out_size,tmp_state.Data(),ps ));
+ boost::shared_ptr<RealSpatialImageState> out_state(new RealSpatialImageState(out_size,ps ));
out_state->SetSpatialOrigin(in_state.GetSpatialOrigin());
out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
@@ -220,21 +164,14 @@ ImageStateBasePtr FFTFnc::VisitState<Complex,HalfFrequencyDomain>(const ComplexH
assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
- OST_FFTW_fftw_complex* fftw_in =
-reinterpret_cast<OST_FFTW_fftw_complex*>(out_ptr);
+ Complex* in_ptr = tmp_state.Data().GetData();
+ OST_FFTW_fftw_complex* fftw_in =reinterpret_cast<OST_FFTW_fftw_complex*>(in_ptr);
- OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft_c2r(rank,n,
- fftw_in,fftw_out,
- FFTW_ESTIMATE);
+ OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft_c2r(rank,n,fftw_in,fftw_out,FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
OST_FFTW_fftw_execute(plan);
OST_FFTW_fftw_destroy_plan(plan);
- // un-pad (leftover from in-place transform)
- Real* dout_ptr = reinterpret_cast<Real*>(out_ptr);
- for(size_t i=1;i<block_count;++i) {
- std::copy(&dout_ptr[i*2*src_size],&dout_ptr[i*2*src_size+dst_size],&dout_ptr[i*dst_size]);
- }
Real fac = 1.0/static_cast<Real>(out_state->GetSize().GetVolume());
for(Real* ptr = out_state->Data().GetData(); ptr<out_state->Data().GetEnd(); ++ptr) {
@@ -249,29 +186,24 @@ reinterpret_cast<OST_FFTW_fftw_complex*>(out_ptr);
template <>
ImageStateBasePtr FFTFnc::VisitState<Complex,SpatialDomain>(const ComplexSpatialImageState& in_state) const
{
+ assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
+
Size size=in_state.GetExtent().GetSize();
PixelSampling ps=in_state.GetSampling();
ps.SetDomain(FREQUENCY);
- boost::shared_ptr<ComplexFrequencyImageState> out_state(new ComplexFrequencyImageState(size,ps));
- out_state->SetSpatialOrigin(in_state.GetSpatialOrigin());
- out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
-
- out_state->Data()=in_state.Data(); // use assignement op to copy data area to new state
int rank = size.GetDim();
int n[3] = {static_cast<int>(size[0]),static_cast<int>(size[1]),static_cast<int>(size[2])};
int dir = FFTW_FORWARD;
+ OST_FFTW_fftw_complex* fftw_in = const_cast<OST_FFTW_fftw_complex*>(reinterpret_cast<const OST_FFTW_fftw_complex*>(in_state.Data().GetData()));
- assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
- OST_FFTW_fftw_complex* fftw_out =
-reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
+ boost::shared_ptr<ComplexFrequencyImageState> out_state(new ComplexFrequencyImageState(size,ps));
+ out_state->SetSpatialOrigin(in_state.GetSpatialOrigin());
+ out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
+ OST_FFTW_fftw_complex* fftw_out = reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
- // in place transform
+ // out of place transform
OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
- OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft(rank,n,
- fftw_out, fftw_out,
- dir,
- FFTW_ESTIMATE);
-
+ OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft(rank,n,fftw_in, fftw_out,dir,FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
OST_FFTW_fftw_execute(plan);
OST_FFTW_fftw_destroy_plan(plan);
@@ -284,35 +216,39 @@ reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
template <>
ImageStateBasePtr FFTFnc::VisitState<Complex,FrequencyDomain>(const ComplexFrequencyImageState& in_state) const
{
- // copy in state
- ComplexFrequencyImageState tmp(in_state);
- if(ori_flag_) tmp.AdjustPhaseOrigin(-in_state.GetSpatialOrigin());
+ assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
Size size=in_state.GetExtent().GetSize();
PixelSampling ps=in_state.GetSampling();
ps.SetDomain(SPATIAL);
- // use memory located for tmp
- boost::shared_ptr<ComplexSpatialImageState> out_state(new ComplexSpatialImageState(size,tmp.Data(),ps));
- out_state->SetSpatialOrigin(in_state.GetSpatialOrigin());
- out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
-
int rank = size.GetDim();
int n[3] = {static_cast<int>(size[0]),static_cast<int>(size[1]),static_cast<int>(size[2])};
int dir = FFTW_BACKWARD;
- assert(sizeof(OST_FFTW_fftw_complex)==sizeof(Complex));
- OST_FFTW_fftw_complex* fftw_out =
-reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
-
- // in place transform
- OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
- OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft(rank,n,
- fftw_out, fftw_out,
- dir,
- FFTW_ESTIMATE);
-
- OST_FFTW_fftw_execute(plan);
- OST_FFTW_fftw_destroy_plan(plan);
+ boost::shared_ptr<ComplexSpatialImageState> out_state(new ComplexSpatialImageState(size,ps));
+ out_state->SetSpatialOrigin(in_state.GetSpatialOrigin());
+ out_state->SetAbsoluteOrigin(in_state.GetAbsoluteOrigin());
+ OST_FFTW_fftw_complex* fftw_out = reinterpret_cast<OST_FFTW_fftw_complex*>(out_state->Data().GetData());
+
+ if(ori_flag_){
+ // copy in state
+ ComplexFrequencyImageState tmp_state(in_state);
+ tmp_state.AdjustPhaseOrigin(-in_state.GetSpatialOrigin());
+ OST_FFTW_fftw_complex* fftw_in = reinterpret_cast<OST_FFTW_fftw_complex*>(tmp_state.Data().GetData());
+ // out of place transform
+ OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
+ OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft(rank,n, fftw_in, fftw_out,dir,FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
+ OST_FFTW_fftw_execute(plan);
+ OST_FFTW_fftw_destroy_plan(plan);
+ }else{
+ // use in state
+ OST_FFTW_fftw_complex* fftw_in = const_cast<OST_FFTW_fftw_complex*>(reinterpret_cast<const OST_FFTW_fftw_complex*>(in_state.Data().GetData()));
+ // out of place transform
+ OST_FFTW_fftw_plan_with_nthreads(std::max<int>(1,IDEAL_NUMBER_OF_THREADS()));
+ OST_FFTW_fftw_plan plan = OST_FFTW_fftw_plan_dft(rank,n, fftw_in, fftw_out,dir,FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
+ OST_FFTW_fftw_execute(plan);
+ OST_FFTW_fftw_destroy_plan(plan);
+ }
Real fac = 1.0/static_cast<Real>(size.GetVolume());
for(Complex* ptr = out_state->Data().GetData(); ptr<out_state->Data().GetEnd(); ++ptr) {
diff --git a/modules/img/alg/src/fractional_shift.hh b/modules/img/alg/src/fractional_shift.hh
index 0328e6620..9d739f302 100644
--- a/modules/img/alg/src/fractional_shift.hh
+++ b/modules/img/alg/src/fractional_shift.hh
@@ -29,7 +29,7 @@
namespace ost { namespace img { namespace alg {
-class DLLEXPORT_IMG_ALG FractionalShift: public ModIPAlgorithm
+class FractionalShift: public ModIPAlgorithm
{
public:
FractionalShift(Real sx=0.0, Real sy=0.0, Real sz=0.0);
diff --git a/modules/img/alg/src/histogram.cc b/modules/img/alg/src/histogram.cc
index 0a094bd83..6225fd99c 100644
--- a/modules/img/alg/src/histogram.cc
+++ b/modules/img/alg/src/histogram.cc
@@ -58,7 +58,7 @@ template <typename T, class D>
void HistogramBase::VisitState(const ImageStateImpl<T,D>& isi)
{
bins_=Bins(bin_count_,0);
- for(T* ptr = isi.Data().GetData(); ptr<isi.Data().GetEnd(); ++ptr) {
+ for(const T* ptr = isi.Data().GetData(); ptr<isi.Data().GetEnd(); ++ptr) {
Real val=Val2Val<T,Real>(*ptr);
val = std::max(min_,val);
val=std::min(max_,val);
diff --git a/modules/img/alg/src/transcendentals.hh b/modules/img/alg/src/transcendentals.hh
index 84de4bd38..60c1f6bcf 100644
--- a/modules/img/alg/src/transcendentals.hh
+++ b/modules/img/alg/src/transcendentals.hh
@@ -36,22 +36,22 @@
#include <ost/img/alg/module_config.hh>
-#define IMG_ALG_TRANSCENDENTALS_BLOCK(FF,NN,SS) \
-struct FF { \
+#define IMG_ALG_TRANSCENDENTALS_BLOCK(FF,NN,SS) \
+struct FF { \
FF() {} \
~FF() {} \
template <typename T, class D> \
void VisitState(ImageStateImpl<T,D>& isi) const { \
- T* end = isi.Data().GetEnd(); \
+ const T* end = isi.Data().GetEnd(); \
for(T* it = isi.Data().GetData(); it!=end; ++it) { \
(*it) = SS (*it); \
} \
} \
template <class D> \
void VisitState(ImageStateImpl<Word,D>& isi) const { \
- Word* end = isi.Data().GetEnd(); \
+ const Word* end = isi.Data().GetEnd(); \
for(Word* it = isi.Data().GetData(); it!=end; ++it) { \
- (*it) = static_cast<Word>(SS(static_cast<Real>(*it))); \
+ (*it) = static_cast<Word>(SS(static_cast<Real>(*it))); \
} \
} \
static String GetAlgorithmName() {return "";} \
@@ -75,7 +75,7 @@ struct PowFnc {
template <typename T, class D>
void VisitState(ImageStateImpl<T,D>& isi) const {
- T* end = isi.Data().GetEnd();
+ const T* end = isi.Data().GetEnd();
for(T* it = isi.Data().GetData(); it!=end; ++it) {
(*it) = std::pow(*it,exp_);
}
diff --git a/modules/img/alg/tests/test_fft.cc b/modules/img/alg/tests/test_fft.cc
index 3e3702690..7d8d0a1d1 100644
--- a/modules/img/alg/tests/test_fft.cc
+++ b/modules/img/alg/tests/test_fft.cc
@@ -408,7 +408,7 @@ void Test_DFT(DataType TYPE)
Real absdiff=std::abs(ri3.GetReal(it)-ri1.GetReal(it));
msg.str("");
msg << "@" << Point(it) << ": abs(" << ri3.GetReal(it) << "-" <<ri1.GetReal(it) << ")="<<absdiff;
- BOOST_REQUIRE(absdiff<1e-5);
+ BOOST_REQUIRE_MESSAGE(absdiff<1e-5,msg.str());
}
}
diff --git a/modules/img/base/src/image_state/image_state_base.hh b/modules/img/base/src/image_state/image_state_base.hh
index 016fab63a..22b5ed14a 100644
--- a/modules/img/base/src/image_state/image_state_base.hh
+++ b/modules/img/base/src/image_state/image_state_base.hh
@@ -57,7 +57,7 @@ public:
// provide deep copy
virtual ImageStateBasePtr Clone(bool cc=true) const = 0;
- virtual long MemSize() const = 0;
+ virtual size_t MemSize() const = 0;
///////////////////////////////////////
// abstract virtual methods
diff --git a/modules/img/base/src/image_state/image_state_impl.cc b/modules/img/base/src/image_state/image_state_impl.cc
index fbef8e2e6..7dec67093 100644
--- a/modules/img/base/src/image_state/image_state_impl.cc
+++ b/modules/img/base/src/image_state/image_state_impl.cc
@@ -117,7 +117,7 @@ ImageStateBasePtr ImageStateImpl<T,D>::Clone(bool cc) const
template <typename T, class D>
-long ImageStateImpl<T,D>::MemSize() const
+size_t ImageStateImpl<T,D>::MemSize() const
{
return data_.MemSize();
}
@@ -125,7 +125,7 @@ long ImageStateImpl<T,D>::MemSize() const
template <typename T, class D>
DataType ImageStateImpl<T,D>::GetType() const
{
- return Val2Type<T>();
+ return data_.GetDataType();
}
template <typename T, class D>
@@ -564,15 +564,15 @@ const T& ImageStateImpl<T,D>::Value(const Index& i) const
}
template <typename T, class D>
-T& ImageStateImpl<T,D>::Value(unsigned int i)
+T& ImageStateImpl<T,D>::Value(size_t i)
{
- return data_.Value(i);
+ return data_.GetData()[i];
}
template <typename T, class D>
-const T& ImageStateImpl<T,D>::Value(unsigned int i) const
+const T& ImageStateImpl<T,D>::Value(size_t i) const
{
- return data_.Value(i);
+ return data_.GetData()[i];
}
template <typename T, class D>
diff --git a/modules/img/base/src/image_state/image_state_impl.hh b/modules/img/base/src/image_state/image_state_impl.hh
index 62e6895cf..151592ee7 100644
--- a/modules/img/base/src/image_state/image_state_impl.hh
+++ b/modules/img/base/src/image_state/image_state_impl.hh
@@ -74,19 +74,6 @@ public:
// also takes absolute center, requires logical extent, not physical one!
ImageStateImpl(const Extent& e, const PixelSampling& s, const Vec3& c);
- /*
- special ctor to utilize pre-allocated memory: the data area allocated
- by ValueHolder will be transfered to this image state
- */
- template <typename U>
- ImageStateImpl(const Extent& e, ValueHolder<U>& d, const PixelSampling& s):
- domain_(e),
- data_(domain_.GetExtent().GetSize(),d),
- sampling_(s)
- {
- assert(d.GetPhysicalSize()*sizeof(U)==data_.GetPhysicalSize()*sizeof(T));
- }
-
virtual ~ImageStateImpl();
SharedPtrType CloneState(bool cc=true) const;
@@ -94,7 +81,7 @@ public:
// image state base interface
virtual ImageStateBasePtr Clone(bool cc=true) const;
- virtual long MemSize() const;
+ virtual size_t MemSize() const;
virtual DataType GetType() const;
@@ -181,10 +168,10 @@ public:
const T& Value(const Index& i) const;
// pass-through to value holder
- T& Value(unsigned int i);
+ T& Value(size_t i);
// pass-through to value holder
- const T& Value(unsigned int i) const;
+ const T& Value(size_t i) const;
//! direct access to value holder
ValueHolder<T>& Data();
diff --git a/modules/img/base/src/image_state/image_state_spatial_domain.hh b/modules/img/base/src/image_state/image_state_spatial_domain.hh
index 5d5e88325..432008e19 100644
--- a/modules/img/base/src/image_state/image_state_spatial_domain.hh
+++ b/modules/img/base/src/image_state/image_state_spatial_domain.hh
@@ -49,8 +49,7 @@ namespace ost { namespace img { namespace image_state {
class DLLEXPORT_OST_IMG_BASE SpatialDomain {
public:
SpatialDomain(const Extent& e):
- extent_(e),
- physical_extent_(ExtractPhysicalExtent(e))
+ extent_(e)
{}
// interface for ImageStateImpl
@@ -73,7 +72,7 @@ public:
}
Extent GetPhysicalExtent() const {
- return physical_extent_;
+ return extent_;
}
template <typename V>
@@ -110,29 +109,13 @@ public:
Index Point2Index(const Point& p) const {
return Index(p[0]-extent_.GetStart()[0],
- p[1]-extent_.GetStart()[1],
- p[2]-extent_.GetStart()[2]);
+ p[1]-extent_.GetStart()[1],
+ p[2]-extent_.GetStart()[2]);
}
private:
Extent extent_;
- Extent physical_extent_;
-
- Extent ExtractPhysicalExtent(const Extent& e) const {
-#if 1
- Point pad;
- if(e.GetDim()==1) {
- pad=Point( (e.GetSize()[0]&0x1) ? 1 : 2,0,0);
- } else if(e.GetDim()==2) {
- pad=Point(0,(e.GetSize()[1]&0x1) ? 1 : 2,0);
- } else {
- pad=Point(0,0,(e.GetSize()[2]&0x1) ? 1 : 2);
- }
- return Extent(e.GetStart(),e.GetEnd()+pad);
-#else
- return Extent(e.GetStart(),e.GetEnd());
-#endif
- }
+
};
}}} // ns
diff --git a/modules/img/base/src/image_state/value_holder.cc b/modules/img/base/src/image_state/value_holder.cc
index ed320e28e..71c79d1d7 100644
--- a/modules/img/base/src/image_state/value_holder.cc
+++ b/modules/img/base/src/image_state/value_holder.cc
@@ -21,7 +21,7 @@
/*
value holder for image states
- Author: Ansgar Philippsen
+ Authors: Ansgar Philippsen, Andreas Schenk
*/
#ifndef IMG_IMAGE_STATE_INST_H
@@ -31,138 +31,98 @@
#include <cassert>
#include <algorithm>
+#include <ost/img/value_util.hh>
#include <ost/img/size.hh>
#include "value_holder.hh"
#include "index.hh"
-#define USE_ROW_ORDER 1
-
namespace ost { namespace img { namespace image_state {
namespace {
-// helper function
-template <class V>
-V* ValueCopy(const V* src, size_t volume_)
-{
- V* dest = new V[volume_]; // if this throws, nothing has happened yet
- try {
- std::copy(src,src+volume_,dest);
- } catch (...) {
- delete []dest; // remove already allocated data
- throw;
- }
- return dest;
-}
-
-template <typename V>
-DataType GetValueType();
-
-template <>
-DataType GetValueType<Complex>() {return COMPLEX;}
-
template <typename V>
-DataType GetValueType() {return REAL;}
-
-}
-
-template <typename V>
-ValueHolder<V>::ValueHolder(unsigned int wi, unsigned int he, unsigned int de):
- width_( (wi<1) ? 1 : wi),
- height_( (he<1) ? 1: he ),
- depth_( (de<1) ? 1: de ),
- volume_(width_*height_*depth_),
- volume2_(volume_),
- data_(0),
- row_(0),
- slice_(0)
+size_t calc_volume(const Size& s)
{
- try {
- data_ = new V[volume2_];
- setup();
- clear();
- } catch(...) {
- delete []data_;
+ unsigned int tmp;
+ unsigned int numbits = 0;
+ //find and add highest bit set for each dimension
+ tmp=s[0];
+ while (tmp >>= 1)
+ {
+ ++numbits;
+ }
+ tmp=s[1];
+ while (tmp >>= 1)
+ {
+ ++numbits;
+ }
+ tmp=s[2];
+ while (tmp >>= 1)
+ {
+ ++numbits;
+ }
+ // check if number of pixels (2**(numbits+1)) larger than what size_t can handle
+ if(numbits>=sizeof(size_t)*sizeof(V)*8){
+ throw std::bad_alloc();
}
+ return static_cast<size_t>(s[0])*static_cast<size_t>(s[1])*static_cast<size_t>(s[2]);
}
+} //ns
+
template <typename V>
ValueHolder<V>::ValueHolder(const Size& s):
- width_(s[0]),
- height_(s[1]),
+#ifdef USE_ROW_ORDER
depth_(s[2]),
- volume_(width_*height_*depth_),
- volume2_(volume_),
- data_(0),
- row_(0),
- slice_(0)
+ height_depth_(s[1]*depth_),
+ data_(s[0]*height_depth_)
+#else
+ height_(s[1]),
+ width_height_(s[0]*height_),
+ data_(calc_volume<V>(s))
+#endif
{
- try {
- data_ = new V[volume2_];
- setup();
- clear();
- } catch(...) {
- delete []data_;
- }
}
template <typename V>
ValueHolder<V>::ValueHolder(const Size& s, const Size& ps):
- width_(s[0]),
- height_(s[1]),
+#ifdef USE_ROW_ORDER
depth_(s[2]),
- volume_(width_*height_*depth_),
- volume2_(ps.GetVolume()),
- data_(0),
- row_(0),
- slice_(0)
+ height_depth_(s[1]*depth_),
+#else
+ height_(s[1]),
+ width_height_(s[0]*height_),
+#endif
+ data_(calc_volume<V>(ps))
{
- try {
- data_ = new V[volume2_];
- setup();
- clear();
- } catch(...) {
- delete []data_;
- }
}
template <typename V>
ValueHolder<V>::ValueHolder(const ValueHolder<V>& h):
- width_(h.width_),
- height_(h.height_),
+#ifdef USE_ROW_ORDER
depth_(h.depth_),
- volume_(h.volume_),
- volume2_(h.volume2_),
- data_(ValueCopy(h.data_,h.volume2_)), // exception safe call
- row_(0),
- slice_(0)
+ height_depth_(h.height_depth_),
+#else
+ height_(h.height_),
+ width_height_(h.width_height_),
+#endif
+ data_(h.data_.begin(),h.data_.end())
{
- try {
- setup();
- } catch(...) {
- delete []data_;
- }
}
template <typename V>
ValueHolder<V>& ValueHolder<V>::operator=(const ValueHolder<V>& h)
{
if(this!=&h) {
- delete []data_;
- delete []row_;
- delete []slice_;
- data_ = ValueCopy(h.data_, h.volume2_);
- width_=h.width_;
- height_=h.height_;
- depth_=h.depth_;
- volume_=h.volume_;
- volume2_=h.volume2_;
- try {
- setup();
- } catch (...) {
- throw;
- }
+ data_.assign(h.data_.begin(),h.data_.end());
+ #ifdef USE_ROW_ORDER
+ depth_=h.depth_;
+ height_depth_=h.height_depth_;
+ #else
+ height_=h.height_;
+ width_height_=h.width_height_;
+ #endif
}
return *this;
}
@@ -170,134 +130,59 @@ ValueHolder<V>& ValueHolder<V>::operator=(const ValueHolder<V>& h)
template <typename V>
ValueHolder<V>::~ValueHolder()
{
- delete []slice_;
- delete []row_;
- delete []data_;
}
template <typename V>
void ValueHolder<V>::Swap(ValueHolder& vh)
{
- std::swap<>(width_,vh.width_);
- std::swap<>(height_,vh.height_);
- std::swap<>(depth_,vh.depth_);
- std::swap<>(volume_,vh.volume_);
- std::swap<>(volume2_,vh.volume2_);
+ #ifdef USE_ROW_ORDER
+ std::swap<>(depth_,vh.depth_);
+ std::swap<>(height_depth_,vh.height_depth_);
+ #else
+ std::swap<>(height_,vh.height_);
+ std::swap<>(width_height_,vh.width_height_);
+ #endif
std::swap<>(data_,vh.data_);
- std::swap<>(row_,vh.row_);
- std::swap<>(slice_,vh.slice_);
}
-template <typename V>
-V* ValueHolder<V>::ReleaseData()
-{
- V* ret = data_; data_=0;
- delete []slice_; slice_=0;
- delete []row_; row_=0;
- return ret;
-}
-
-template <typename V>
-Size ValueHolder<V>::GetSize() const
-{
- return Size(width_,depth_,height_);
-}
template <typename V>
V& ValueHolder<V>::Value(const Index& i)
{
- assert(i.u<width_ && i.v<height_ && i.w<depth_);
#ifdef USE_ROW_ORDER
- //return ((slice_[i.u])[i.v])[i.w];
- return data_[i.u*height_*depth_ + i.v*depth_ +i.w];
+ assert(i.w<depth_);
+ return data_[i.u*height_depth_ + i.v*depth_ +i.w];
#else
- //return ((slice_[i.w])[i.v])[i.u];
- return data_[i.w*width_*height_ + i.v*height_ +i.u];
+ assert(i.v<height_);
+ return data_[i.w*width_height_ + i.v*height_ +i.u];
#endif
}
template <typename V>
const V& ValueHolder<V>::Value(const Index& i) const
{
- assert(i.u<width_ && i.v<height_ && i.w<depth_);
#ifdef USE_ROW_ORDER
- //return ((slice_[i.u])[i.v])[i.w];
- return data_[i.u*height_*depth_ + i.v*depth_ +i.w];
+ assert(i.w<depth_);
+ return data_[i.u*height_depth_ + i.v*depth_ +i.w];
#else
- //return ((slice_[i.w])[i.v])[i.u];
- return data_[i.w*width_*height_ + i.v*height_ +i.u];
+ assert(i.v<height_);
+ return data_[i.w*width_height_ + i.v*height_ +i.u];
#endif
}
-template <typename V>
-V& ValueHolder<V>::Value(unsigned int i)
-{
- assert(i<volume_);
- return data_[i];
-}
-
-template <typename V>
-const V& ValueHolder<V>::Value(unsigned int i) const
-{
- assert(i<volume_);
- return data_[i];
-}
template <typename V>
DataType ValueHolder<V>::GetDataType()
{
- return GetValueType<V>();
+ return Val2Type<V>();
}
template <typename V>
-long ValueHolder<V>::MemSize() const
+size_t ValueHolder<V>::MemSize() const
{
- //return sizeof(V)*volume2_+sizeof(V*)*width_*height_+sizeof(V**)*width_+sizeof(*this);
- return sizeof(V)*volume2_;
+ return sizeof(V)*data_.size();
}
-template <typename V>
-void ValueHolder<V>::setup()
-{
- typedef V* VPtr;
- typedef VPtr* VPtrPtr;
- row_=0;
- slice_=0;
- try {
-#if 0
-#ifdef USE_ROW_ORDER
- row_ = new VPtr[width_*height_];
- slice_ = new VPtrPtr[width_];
- for(unsigned int u=0;u<width_;++u) {
- for(unsigned int v=0;v<height_;++v) {
- row_[u*height_+v]=&data_[u*height_*depth_+v*depth_];
- }
- slice_[u]=&row_[u*height_];
- }
-#else
- row_ = new VPtr[height_*depth_];
- slice_ = new VPtrPtr[depth_];
-
- for(unsigned int w=0;w<depth_;++w) {
- for(unsigned int v=0;v<height_;++v) {
- row_[w*height_+v]=&data_[w*width_*height_+v*width_];
- }
- slice_[w]=&row_[w*height_];
- }
-#endif
-#endif
- } catch (...) {
- delete slice_;
- delete row_;
- throw;
- }
-}
-
-template <typename V>
-void ValueHolder<V>::clear()
-{
- for(unsigned int c=0;c<volume2_;++c) data_[c]=V();
-}
}}} // namespaces
diff --git a/modules/img/base/src/image_state/value_holder.hh b/modules/img/base/src/image_state/value_holder.hh
index e0692f20a..6caa155cc 100644
--- a/modules/img/base/src/image_state/value_holder.hh
+++ b/modules/img/base/src/image_state/value_holder.hh
@@ -21,7 +21,7 @@
/*
value holder for image state
- Author: Ansgar Philippsen
+ Authors: Ansgar Philippsen, Andreas Schenk
*/
#ifndef VALUE_HOLDER_H
@@ -41,6 +41,10 @@ void ReleaseAndReconstruct();
}
+#define USE_ROW_ORDER 1
+
+
+
namespace ost { namespace img { namespace image_state {
/*
@@ -57,14 +61,12 @@ template <typename V>
class TEMPLATE_EXPORT ValueHolder {
public:
typedef V* VPtr;
+ typedef const V* ConstVPtr;
typedef VPtr* VPtrPtr;
/*! @name Construction, Release, etc
*/
//@{
- //! initialization with explicit width, height, and depth
- ValueHolder(unsigned int wi, unsigned int he, unsigned int de);
-
//! initialization with size
ValueHolder(const Size& s);
@@ -74,25 +76,6 @@ public:
//! copy ctor provides full copy!
ValueHolder(const ValueHolder<V>& h);
- //! initialize with a size and an existing value holder, grabbing the memory
- template <typename U>
- ValueHolder(const Size& s, ValueHolder<U>& v):
- width_(s[0]),
- height_(s[1]),
- depth_(s[2]),
- volume_(width_*height_*depth_),
- volume2_(v.GetPhysicalSize()*sizeof(U)/sizeof(V)),
- data_(reinterpret_cast<V*>(v.ReleaseData())),
- row_(0),
- slice_(0)
- {
- try {
- setup();
- } catch(...) {
- delete []data_;
- }
- }
-
//! assignement provides full copy!
/*!
@@ -105,125 +88,52 @@ public:
//! free allocated memory upon destruction
~ValueHolder();
- //! release the held data, invalidates holder!
- /*!
- the data area is returned, but no longer belongs
- to the value holder. it now belongs to the caller, who
- must eiter free it explicitely, or use it in a
- constructor argument to a new value holder.
- */
- VPtr ReleaseData();
//! swap data with another value holder
void Swap(ValueHolder& vh);
//@}
- //! @name Property access
- //@{
-
- // retrieve width, at least 1
- unsigned int GetWidth() const {return width_;}
- // retrieve height, at least 1
- unsigned int GetHeight() const {return height_;}
- // retrieve depth, at least 1
- unsigned int GetDepth() const {return depth_;}
- // retrieve overall size (width*height*depth), at least 1
- unsigned int GetVolume() const {return volume_;}
- // retrieve dimensions as Size object
- Size GetSize() const;
-
static DataType GetDataType();
- unsigned int GetPhysicalSize() const {return volume2_;}
-
- long MemSize() const;
+ size_t MemSize() const;
- //@}
- /*! @name Data Access
- access is either via interface that uses the internal
- row and slice pointers, or these pointers may be
- retrieved directly.
- */
//! return direct r/w access to the value without boundary check
/*!
The lookup is based on an integer triplet encapsulated within
- Index. Efficient lookup using helper tables.
+ Index.
*/
V& Value(const Index& i);
//! return direct ro access to the value without boundary check
/*!
The lookup is based on an integer triplet encapsulated within
- Index. Efficient lookup using helper tables.
+ Index.
*/
const V& Value(const Index& i) const;
- //! return direct r/w access to the value without boundary check
- /*!
- DEPRECATED! Use GetData()[i]
- */
- V& Value(unsigned int i);
-
- //! return direct ro access to the value without boundary check
- /*!
- DEPRECATED! Use GetData()[i]
- */
- const V& Value(unsigned int i) const;
-
-
//! return pointer to raw data
- VPtr GetData() {return data_;}
+ VPtr GetData() {return &data_[0];}
//! const version of GetData()
- const VPtr GetData() const {return data_;}
-
- //! return number of data items DEPRECATED
- int DataCount() const {return volume_;}
-
- //! return number of data items
- /*!
- this is the volume: width*height*depth
- */
- int GetDataCount() const {return volume_;}
+ ConstVPtr GetData() const {return &data_[0];}
- const VPtr GetEnd() const {return &data_[volume_];}
- //! return pointer to row pointers
- VPtr* GetRows() {return row_;}
- //! const version of GetRows()
- const VPtr* GetRows() const {return row_;}
- //! return number of row pointers
- /*!
- this is width*height
- */
- int GetRowCount() const {return width_*height_;}
+ ConstVPtr GetEnd() const {return &data_[0]+data_.size();}
- //! return pointer to slice pointers
- VPtrPtr* GetSlices() {return slice_;}
- //! const version of GetSlices()
- const VPtrPtr* GetSlices() const {return slice_;}
- //! return number of slices
- /*!
- equals the width
- */
- int GetSliceCount() const {return width_;}
- //@}
private:
- unsigned int width_, height_, depth_, volume_;
- // this is a hack to get padding for fftw to work...
- unsigned int volume2_;
+#ifdef USE_ROW_ORDER
+ size_t depth_;
+ size_t height_depth_;
+#else
+ size_t height_;
+ size_t width_height_;
+#endif
// actual data storage
- V* data_;
- // optimization
- V** row_;
- // optimization
- V*** slice_;
-
- void setup();
- void clear();
+ std::vector<V> data_;
+
};
}}} // namespaces
diff --git a/modules/img/base/tests/test_domains.cc b/modules/img/base/tests/test_domains.cc
index 4fdc4ad25..00357322c 100644
--- a/modules/img/base/tests/test_domains.cc
+++ b/modules/img/base/tests/test_domains.cc
@@ -43,31 +43,31 @@ void test_spatial()
SpatialDomain sd(Extent(Size(5,6)));
Extent le(Point(0,0),Point(4,5));
- Extent pe(Point(0,0),Point(4,7));
+ Extent pe(Point(0,0),Point(4,5));
BOOST_CHECK(sd.GetDomain()==SPATIAL);
mesg.str("");
mesg << " expected " << le << " but got " << sd.GetExtent();
- BOOST_CHECK(sd.GetExtent()==le);
+ BOOST_CHECK_MESSAGE(sd.GetExtent()==le,mesg.str());
mesg.str("");
mesg << " expected " << le << " but got " << sd.GetLogicalExtent();
- BOOST_CHECK(sd.GetLogicalExtent()==le);
+ BOOST_CHECK_MESSAGE(sd.GetLogicalExtent()==le,mesg.str());
mesg.str("");
mesg << " expected " << pe << " but got " << sd.GetPhysicalExtent();
- BOOST_CHECK(sd.GetPhysicalExtent()==pe);
+ BOOST_CHECK_MESSAGE(sd.GetPhysicalExtent()==pe,mesg.str());
sd=SpatialDomain(Extent(Size(3,4,5)));
le=Extent(Point(0,0,0),Point(2,3,4));
- pe=Extent(Point(0,0,0),Point(2,3,5));
+ pe=Extent(Point(0,0,0),Point(2,3,4));
mesg.str("");
mesg << " expected " << le << " but got " << sd.GetExtent();
- BOOST_CHECK(sd.GetExtent()==le);
+ BOOST_CHECK_MESSAGE(sd.GetExtent()==le,mesg.str());
mesg.str("");
mesg << " expected " << le << " but got " << sd.GetLogicalExtent();
- BOOST_CHECK(sd.GetLogicalExtent()==le);
+ BOOST_CHECK_MESSAGE(sd.GetLogicalExtent()==le,mesg.str());
mesg.str("");
mesg << " expected " << pe << " but got " << sd.GetPhysicalExtent();
- BOOST_CHECK(sd.GetPhysicalExtent()==pe);
+ BOOST_CHECK_MESSAGE(sd.GetPhysicalExtent()==pe,mesg.str());
}
void test_half_freq()
@@ -89,11 +89,6 @@ void test_half_freq()
mesg << " expected " << pe << " but got " << hd.GetPhysicalExtent();
BOOST_CHECK_MESSAGE(hd.GetPhysicalExtent()==pe,mesg.str());
- SpatialDomain sd(Extent(Size(5,6)));
- mesg.str("");
- mesg << " expected volume match " << hd.GetPhysicalExtent() << " " << sd.GetPhysicalExtent();
- BOOST_CHECK_MESSAGE(hd.GetPhysicalExtent().GetVolume()*2==sd.GetPhysicalExtent().GetVolume(),mesg.str());
-
// 3d
hd=HalfFrequencyDomain(Extent(Size(3,4,5)));
le=Extent(Point(-1,-1,-2),Point(1,2,2));
@@ -110,10 +105,6 @@ void test_half_freq()
mesg << " expected " << pe << " but got " << hd.GetPhysicalExtent();
BOOST_CHECK_MESSAGE(hd.GetPhysicalExtent()==pe,mesg.str());
- sd=SpatialDomain(Extent(Size(3,4,5)));
- mesg.str("");
- mesg << " expected volume match " << hd.GetPhysicalExtent() << " " << sd.GetPhysicalExtent();
- BOOST_CHECK_MESSAGE(hd.GetPhysicalExtent().GetVolume()*2==sd.GetPhysicalExtent().GetVolume(),mesg.str());
}
diff --git a/modules/img/base/tests/test_value_holder.cc b/modules/img/base/tests/test_value_holder.cc
index d859912dd..f5b09f938 100644
--- a/modules/img/base/tests/test_value_holder.cc
+++ b/modules/img/base/tests/test_value_holder.cc
@@ -35,15 +35,6 @@ using namespace ost::img::image_state;
namespace value_holder_test {
-template <typename V>
-void Construction()
-{
- ValueHolder<V> v(8,6,10);
- BOOST_CHECK (v.GetWidth()==8);
- BOOST_CHECK (v.GetHeight()==6);
- BOOST_CHECK (v.GetDepth()==10);
- BOOST_CHECK (v.GetVolume()==8*6*10);
-}
template <typename V>
void ReadWrite()
@@ -51,7 +42,7 @@ void ReadWrite()
int wi=4;
int he=3;
int de=5;
- ValueHolder<V> vh(wi,he,de);
+ ValueHolder<V> vh(Size(wi,he,de));
for(int u=0;u<wi;++u) {
for(int v=0;v<he;++v) {
@@ -60,7 +51,7 @@ void ReadWrite()
int indx = w+de*(v+he*u);
vh.Value(Index(u,v,w))=val;
BOOST_REQUIRE(vh.Value(Index(u,v,w)) == val);
- BOOST_REQUIRE(vh.Value(indx) == val);
+ BOOST_REQUIRE(vh.GetData()[indx] == val);
}
}
}
@@ -73,7 +64,7 @@ void CopyCtor()
int he=3;
int de=5;
- ValueHolder<V> vh1(wi,he,de);
+ ValueHolder<V> vh1(Size(wi,he,de));
for(int u=0;u<wi;++u) {
for(int v=0;v<he;++v) {
@@ -104,8 +95,6 @@ test_suite* CreateValueHolderTest()
test_suite* ts=BOOST_TEST_SUITE("ValueHolder Test");
- ts->add( BOOST_TEST_CASE(Construction<Real>) );
- ts->add( BOOST_TEST_CASE(Construction<Complex>) );
ts->add( BOOST_TEST_CASE(ReadWrite<Real>) );
ts->add( BOOST_TEST_CASE(ReadWrite<Complex>) );
ts->add( BOOST_TEST_CASE(CopyCtor<Real>) );
diff --git a/modules/io/src/img/map_io_mrc_handler.cc b/modules/io/src/img/map_io_mrc_handler.cc
index b49e976ec..587958bf0 100644
--- a/modules/io/src/img/map_io_mrc_handler.cc
+++ b/modules/io/src/img/map_io_mrc_handler.cc
@@ -806,7 +806,7 @@ void import_helper(img::MapHandle& image, std::istream& in,const MRC& formatmrc)
static_cast<Real>(header.y)/static_cast<Real>(header.ny),
static_cast<Real>(header.z)/static_cast<Real>(header.nz)));
}else{
- LOG_INFO("Suspicious dell dimensions found. Cannot set sampling.");
+ LOG_INFO("Suspicious cell dimensions found. Cannot set sampling.");
}
LOG_INFO("resulting image extent: " << image.GetExtent());
if(img::image_state::RealSpatialImageState *rs=dynamic_cast<img::image_state::RealSpatialImageState*>(image.ImageStatePtr().get())) {
--
GitLab