Skip to content

Commit

Permalink
implemented DataContainer.min(), fixes #1233
Browse files Browse the repository at this point in the history
  • Loading branch information
evgueni-ovtchinnikov committed Feb 20, 2024
1 parent a323b88 commit 8697ea3
Show file tree
Hide file tree
Showing 12 changed files with 175 additions and 0 deletions.
14 changes: 14 additions & 0 deletions src/Registration/cReg/NiftiImageData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1773,6 +1773,20 @@ void NiftiImageData<dataType>::max(void* ptr) const
*ptr_s = s;
}

template<class dataType>
void NiftiImageData<dataType>::min(void* ptr) const
{
unsigned i = 0;
float s = _data[i++];
for (; i < this->_nifti_image->nvox; ++i) {
float si = _data[i];
if (si < s)
s = si;
}
float* ptr_s = static_cast<float*>(ptr);
*ptr_s = s;
}

template<class dataType>
void NiftiImageData<dataType>::axpby(
const void* ptr_a, const DataContainer& a_x,
Expand Down
1 change: 1 addition & 0 deletions src/Registration/cReg/include/sirf/Reg/NiftiImageData.h
Original file line number Diff line number Diff line change
Expand Up @@ -624,6 +624,7 @@ class NiftiImageData : public ImageData
/// below all void* are actually float*
virtual void sum (void* ptr) const;
virtual void max (void* ptr) const;
virtual void min (void* ptr) const;
virtual void dot (const DataContainer& a_x, void* ptr) const;
virtual void axpby (const void* ptr_a, const DataContainer& a_x, const void* ptr_b, const DataContainer& a_y);
virtual void xapyb (const DataContainer& a_x, const void* ptr_a, const DataContainer& a_y, const void* ptr_b);
Expand Down
10 changes: 10 additions & 0 deletions src/common/SIRF.py
Original file line number Diff line number Diff line change
Expand Up @@ -286,6 +286,16 @@ def max(self):
return z[0]
return z[0] + 1j*z[1]

def min(self):
'''
Returns the maximum of the elements of self data
'''
z = numpy.zeros((2,), dtype=numpy.float32)
try_calling(pysirf.cSIRF_compute_min(self.handle, z.ctypes.data))
if z[1] == 0:
return z[0]
return z[0] + 1j*z[1]

def add(self, other, out=None):
'''
Addition for data containers.
Expand Down
4 changes: 4 additions & 0 deletions src/common/Utilities.py
Original file line number Diff line number Diff line change
Expand Up @@ -604,6 +604,10 @@ def test_data_container_algebra(test, x, eps=1e-5):
t = numpy.max(ax)
test.check_if_equal_within_tolerance(t, s, 0, eps);

s = x.min()
t = numpy.min(ax)
test.check_if_equal_within_tolerance(t, s, 0, eps);

s = x.sum()
t = numpy.sum(ax)
r = numpy.sum(abs(ax))
Expand Down
12 changes: 12 additions & 0 deletions src/common/csirf.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -178,6 +178,18 @@ cSIRF_compute_max(const void* ptr_x, void* ptr_z)
CATCH;
}

extern "C"
void*
cSIRF_compute_min(const void* ptr_x, void* ptr_z)
{
try {
auto const& x = objectFromHandle<DataContainer>(ptr_x);
x.min(ptr_z);
return new DataHandle;
}
CATCH;
}

extern "C"
void*
cSIRF_axpby(
Expand Down
3 changes: 3 additions & 0 deletions src/common/include/sirf/common/DataContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,9 @@ namespace sirf {
/// calculates the value of this container's element with the largest real part
virtual void max(void* ptr) const = 0;

/// calculates the value of this container's element with the smallest real part
virtual void min(void* ptr) const = 0;

/// \c *this = the elementwise product \c x*y
virtual void multiply
(const DataContainer& x, const DataContainer& y) = 0;
Expand Down
1 change: 1 addition & 0 deletions src/common/include/sirf/common/csirf.h
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ void* cSIRF_norm(const void* ptr_x);
void* cSIRF_compute_dot(const void* ptr_x, const void* ptr_y, PTR_FLOAT ptr_z);
void* cSIRF_compute_sum(const void* ptr_x, PTR_FLOAT ptr_z);
void* cSIRF_compute_max(const void* ptr_x, PTR_FLOAT ptr_z);
void* cSIRF_compute_min(const void* ptr_x, PTR_FLOAT ptr_z);
void* cSIRF_axpby(const PTR_FLOAT ptr_a, const void* ptr_x,
const PTR_FLOAT ptr_b, const void* ptr_y);
void* cSIRF_axpbyAlt(const PTR_FLOAT ptr_a, const void* ptr_x,
Expand Down
59 changes: 59 additions & 0 deletions src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -579,6 +579,24 @@ MRAcquisitionData::max(const ISMRMRD::Acquisition& acq_a)
return z;
}

complex_float_t
MRAcquisitionData::min(const ISMRMRD::Acquisition& acq_a)
{
const complex_float_t* pa;
complex_float_t z = 0;
bool init = true;
for (pa = acq_a.data_begin(); pa != acq_a.data_end(); pa++) {
complex_float_t zi = *pa;
float r = std::real(z);
float ri = std::real(zi);
if (init || ri < r) {
z = zi;
init = false;
}
}
return z;
}

void
MRAcquisitionData::dot(const DataContainer& dc, void* ptr) const
{
Expand Down Expand Up @@ -648,6 +666,31 @@ MRAcquisitionData::max(void* ptr) const
*ptr_z = z;
}

void
MRAcquisitionData::min(void* ptr) const
{
int n = number();
complex_float_t z = 0;
ISMRMRD::Acquisition a;
bool init = true;
for (int i = 0; i < n;) {
if (!get_acquisition(i, a)) {
i++;
continue;
}
complex_float_t zi = MRAcquisitionData::min(a);
float r = std::real(z);
float ri = std::real(zi);
if (init || ri < r) {
z = zi;
init = false;
}
i++;
}
complex_float_t* ptr_z = static_cast<complex_float_t*>(ptr);
*ptr_z = z;
}

void
MRAcquisitionData::axpby(
const void* ptr_a, const DataContainer& a_x,
Expand Down Expand Up @@ -1472,6 +1515,22 @@ GadgetronImageData::max(void* ptr) const
*ptr_z = z;
}

void
GadgetronImageData::min(void* ptr) const
{
complex_float_t z = 0;
for (unsigned int i = 0; i < number(); i++) {
const ImageWrap& wi = image_wrap(i);
complex_float_t zi = wi.min();
float r = std::real(z);
float ri = std::real(zi);
if (i == 0 || ri < r)
z = zi;
}
complex_float_t* ptr_z = static_cast<complex_float_t*>(ptr);
*ptr_z = z;
}

void
GadgetronImageData::axpby(
const void* ptr_a, const DataContainer& a_x,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,7 @@ namespace sirf {
static complex_float_t sum(const ISMRMRD::Acquisition& acq_x);
// the value of the element of x with the largest real part
static complex_float_t max(const ISMRMRD::Acquisition& acq_x);
static complex_float_t min(const ISMRMRD::Acquisition& acq_x);
// elementwise multiplication
// y := x .* y
static void multiply
Expand Down Expand Up @@ -601,6 +602,7 @@ namespace sirf {
/// below all void* are actually complex_float_t*
virtual void sum(void* ptr) const;
virtual void max(void* ptr) const;
virtual void min(void* ptr) const;
virtual void dot(const DataContainer& dc, void* ptr) const;
complex_float_t dot(const DataContainer& a_x)
{
Expand Down Expand Up @@ -956,6 +958,7 @@ namespace sirf {
/// below all void* are actually complex_float_t*
virtual void sum(void* ptr) const;
virtual void max(void* ptr) const;
virtual void min(void* ptr) const;
virtual void dot(const DataContainer& dc, void* ptr) const;
virtual void axpby(
const void* ptr_a, const DataContainer& a_x,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -557,6 +557,12 @@ namespace sirf {
IMAGE_PROCESSING_SWITCH_CONST(type_, max_, ptr_, &s);
return s;
}
complex_float_t min() const
{
complex_float_t s;
IMAGE_PROCESSING_SWITCH_CONST(type_, min_, ptr_, &s);
return s;
}
float diff(ImageWrap& iw) const
{
float s;
Expand Down Expand Up @@ -975,6 +981,22 @@ namespace sirf {
}
}

template<typename T>
void min_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
{
const T* i;
*s = 0;
size_t ii = 0;
size_t n = ptr->getNumberOfDataElements();
for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
complex_float_t si = (complex_float_t)*i;
float r = std::real(*s);
float ri = std::real(si);
if (ii == 0 || ri < r)
*s = si;
}
}

template<typename T>
void diff_(const ISMRMRD::Image<T>* ptr_im, float *s) const
{
Expand Down
2 changes: 2 additions & 0 deletions src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,7 @@ namespace sirf {
/// below all void* are actually float*
virtual void sum(void* ptr) const;
virtual void max(void* ptr) const;
virtual void min(void* ptr) const;
virtual void dot(const DataContainer& a_x, void* ptr) const;
float dot(const DataContainer& a_x) const
{
Expand Down Expand Up @@ -1072,6 +1073,7 @@ namespace sirf {
/// below all void* are actually float*
virtual void sum(void* ptr) const;
virtual void max(void* ptr) const;
virtual void min(void* ptr) const;
virtual void dot(const DataContainer& a_x, void* ptr) const;
virtual void axpby(
const void* ptr_a, const DataContainer& a_x,
Expand Down
44 changes: 44 additions & 0 deletions src/xSTIR/cSTIR/stir_data_containers.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,34 @@ STIRAcquisitionData::max(void* ptr) const
*ptr_t = (float)t;
}

void
STIRAcquisitionData::min(void* ptr) const
{
int n = get_max_segment_num();
float t = 0;
bool init = true;
TOF_LOOP
for (int s = 0; s <= n; ++s)
{
SegmentBySinogram<float> seg = get_segment_by_sinogram(s TOF_ARG);
SegmentBySinogram<float>::full_iterator seg_iter;
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
if (init) {
t = *seg_iter++;
init = false;
}
else
t = std::min(t, *seg_iter++);
if (s != 0) {
seg = get_segment_by_sinogram(-s TOF_ARG);
for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
t = std::min(t, *seg_iter++);
}
}
float* ptr_t = static_cast<float*>(ptr);
*ptr_t = (float)t;
}

void
STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const
{
Expand Down Expand Up @@ -496,6 +524,22 @@ STIRImageData::max(void* ptr) const
*ptr_s = (float)s;
}

void
STIRImageData::min(void* ptr) const
{
#if defined(_MSC_VER) && _MSC_VER < 1900
Image3DF::const_full_iterator iter;
#else
typename Array<3, float>::const_full_iterator iter;
#endif
iter = data().begin_all();
float s = *iter++;
for (; iter != data().end_all(); iter++)
s = std::min(s, *iter);
float* ptr_s = static_cast<float*>(ptr);
*ptr_s = (float)s;
}

void
STIRImageData::dot(const DataContainer& a_x, void* ptr) const
{
Expand Down

0 comments on commit 8697ea3

Please sign in to comment.