diff --git a/modules/img/alg/src/fft.cc b/modules/img/alg/src/fft.cc
index ddc137d1a5272ee98d091b211458ce3791f0c7c7..955f218346844dedbb6e816606664be9f864169b 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 0328e66202234d16f8d2cca4d85b54a91696bf58..9d739f3021623ddab8b25e2bfbc6fa5c54ea59e0 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 0a094bd8319675c59f587a7ea1638c6f79ed78c1..6225fd99c017de03e6795c4232eb669b0896b9cb 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 84de4bd3812e4c7996bc5c75ad0b0f061617fd8d..60c1f6bcfd24bc7f759090c4c62d48c926745308 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 3e3702690117843ed71fc633b0ddbfeb69e0d418..7d8d0a1d125b91424af099b77db3c6f1a53d3a13 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 016fab63a063394b5515e42dbb482bf302b874cc..22b5ed14ae02d04049f129456d417a012c3a285e 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 fbef8e2e68dcb65f89f3f8dcd74f9176168da83f..7dec67093909b1536d96a72051084581793969bd 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 62e6895cf03458c64d5798c39017a992271e0211..151592ee781fecfba9f87b83c4836137ea67a5f1 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 5d5e883252984519a3ca5d52ebf348a7000bcd7a..432008e194c3e9bf045bcab2e3fcef5043ceac4a 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 ed320e28e91c44a1232e25c1cd96aa687511ee06..71c79d1d73cd384684669b9f349e1f3b41443d6e 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 e0692f20af8a7810d626fc18aef40b491b134b3d..6caa155cc9fe254dd962d198685347aa60de1e22 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 4fdc4ad2548dd366e32b6598204b927413f5fc93..00357322ca3141f212d7e94f261c32235098527f 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 d859912dd471c774ab21ef842e8a013b350eee5d..f5b09f938f0bd48ba09e0701631d10ab67f4e8c4 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 b49e976eca8e9d1b45aa24b8cffa347934fb8dfe..587958bf0f0c8ea99ade7da9037c1f5addca7ffc 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())) {