Skip to content

Commit

Permalink
Merge pull request #39 from OpenPIV/34-investigate-alternate-fft-impl…
Browse files Browse the repository at this point in the history 
…ementations

modify real cross correlate, add timing for processing
  • Loading branch information
@timdewhirst
timdewhirst authored Nov 28, 2022
2 parents 50b9b84 + 54939b0 commit bd58d21
Show file tree
Hide file tree
Showing 6 changed files with 150 additions and 50 deletions.
19 changes: 14 additions & 5 deletions examples/process/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,28 +150,28 @@ int main( int argc, char* argv[] )
std::vector<point_vector> found_peaks( grid.size() );

// wrap correlators
using correlator_t = std::function<core::cf_image(const core::gf_image&, const core::gf_image&)>;
using correlator_t = std::function<core::gf_image(const core::gf_image&, const core::gf_image&)>;
std::unordered_map<std::string, correlator_t> correlators = {
{"complex",
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::cf_image
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::gf_image
{
static algos::FFT fft{ ia };
return fft.cross_correlate(im_a, im_b);
} },
{"real",
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::cf_image
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::gf_image
{
static algos::FFT fft{ ia };
return fft.cross_correlate_real(im_a, im_b);
} },
{"pocket",
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::cf_image
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::gf_image
{
static algos::PocketFFT fft{ ia };
return fft.cross_correlate(im_a, im_b);
} },
{"pocket_real",
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::cf_image
[ia](const core::gf_image& im_a, const core::gf_image& im_b) -> core::gf_image
{
static algos::PocketFFT fft{ ia };
return fft.cross_correlate_real(im_a, im_b);
Expand Down Expand Up @@ -236,6 +236,8 @@ int main( int argc, char* argv[] )
found_peaks[i] = std::move(result);
};

const auto t1 = std::chrono::high_resolution_clock::now();

// check execution
if (thread_count <= 1)
{
Expand Down Expand Up @@ -299,6 +301,13 @@ int main( int argc, char* argv[] )
}
}

const auto t2 = std::chrono::high_resolution_clock::now();
const std::chrono::duration<double, std::micro> total_us = t2 - t1;
logger::info(
"processing time: {}us, {}us per interrogation area",
total_us,
total_us/found_peaks.size());

// dump output
for ( const auto& pv : found_peaks )
std::cout << pv.xy[0] << ", " << pv.xy[1] << ", " << pv.vxy[0] << ", " << pv.vxy[1] << ", " << pv.sn << "\n";
Expand Down
26 changes: 16 additions & 10 deletions openpiv/algos/fft.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,37 +187,43 @@ namespace openpiv::algos {

template < template <typename> class ImageT,
typename ContainedT,
typename ValueT = typename ContainedT::value_t,
typename OutT = image<g<ValueT>>,
typename = typename std::enable_if_t< is_imagetype_v<ImageT<ContainedT>> >
>
const cf_image& cross_correlate( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
OutT
cross_correlate( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
{
cf_image a_fft{ transform( a, direction::FORWARD ) };
const cf_image& b_fft = transform( b, direction::FORWARD );

a_fft = b_fft * conj( a_fft );
cache().output = real( transform( a_fft, direction::REVERSE ) );
swap_quadrants( cache().output );
OutT output{ real( transform( a_fft, direction::REVERSE ) ) };
swap_quadrants( output );

return cache().output;
return output;
}

template < template <typename> class ImageT,
typename ContainedT,
typename ValueT = typename ContainedT::value_t,
typename OutT = image<g<ValueT>>,
typename = typename std::enable_if_t<
is_imagetype_v<ImageT<ContainedT>> &&
is_real_mono_pixeltype_v<ContainedT>
>
>
const cf_image& cross_correlate_real( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
OutT
cross_correlate_real( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
{
auto [a_fft, b_fft] = transform_real( a, b, direction::FORWARD );
a_fft = b_fft * conj( a_fft );
cache().output = real( transform( a_fft, direction::REVERSE ) );
swap_quadrants( cache().output );
OutT output{ real( transform( a_fft, direction::REVERSE ) ) };
swap_quadrants( output );

return cache().output;
return output;
}

template < template <typename> class ImageT,
Expand Down
133 changes: 99 additions & 34 deletions openpiv/algos/pocket_fft.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ namespace openpiv::algos {
is_real_mono_pixeltype_v<ContainedT>
>
>
std::tuple<cf_image, cf_image>
std::tuple<cf_image&, cf_image&>
transform_real( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b,
direction d = direction::FORWARD ) const
Expand All @@ -144,69 +144,134 @@ namespace openpiv::algos {
<< ", " << size_;
}

// copy data to (real, imag), converting to complex
cache().output = transform( join_from_channels(a, b), d );
using value_t = typename ContainedT::value_t;

// and unravel
const auto& transformed = cache().output;
auto out_a = cf_image( transformed.size() );
auto out_b = cf_image( transformed.size() );
cache().output.resize( a.size() );
cache().temp.resize( b.size() );

const auto width = transformed.width();
const auto height = transformed.height();
for ( uint32_t h=1; h<height/2; ++h)
{
for (uint32_t w=1; w<width; ++w)
auto& out_a = cache().output;
auto& out_b = cache().temp;

constexpr auto stride_lambda = [](auto& im) -> pfft::stride_t
{
const auto t1 = transformed[ {w, h} ];
const auto t2 = transformed[ {width - w, height - h} ];
auto a = 0.5*(t1 + t2.conj());
out_a[ {w, h} ] = a;
out_a[ {width - w, height - h} ] = a.conj();

auto b = 0.5*(t1 - t2.conj());
b = c_f(b.imag, -1.0 * b.real);
out_b[ {w, h} ] = b;
out_b[ {width - w, height - h} ] = b.conj();
}
}
const auto [stride_x, stride_y] = im.stride();
return {static_cast<long>(stride_x), static_cast<long>(stride_y)};
};

const pfft::shape_t shape = {size_.width(), size_.height()};
const pfft::stride_t in_a_stride = stride_lambda(a);
const pfft::stride_t in_b_stride = stride_lambda(b);
const pfft::stride_t out_a_stride = stride_lambda(out_a);
const pfft::stride_t out_b_stride = stride_lambda(out_b);

pfft::r2c<value_t>(
shape,
in_a_stride,
out_a_stride,
{ 0, 1 }, // axes
d == direction::FORWARD, // forward
reinterpret_cast<const value_t*>(a.data()),
reinterpret_cast<std::complex<value_t>*>(out_a.data()),
1.0 );

pfft::r2c<value_t>(
shape,
in_b_stride,
out_b_stride,
{ 0, 1 }, // axes
d == direction::FORWARD, // forward
reinterpret_cast<const value_t*>(b.data()),
reinterpret_cast<std::complex<value_t>*>(out_b.data()),
1.0 );

return { out_a, out_b };
}

template < template <typename> class ImageT,
typename ContainedT,
typename ValueT = typename ContainedT::value_t,
typename OutT = image<g<ValueT>>,
typename = typename std::enable_if_t<
is_imagetype_v<ImageT<ContainedT>> &&
is_complex_mono_pixeltype_v<ContainedT>
>
>
OutT
transform_real( const ImageT<ContainedT>& in,
direction d = direction::FORWARD ) const
{
DECLARE_ENTRY_EXIT
if ( in.size() != size_ )
{
exception_builder< std::runtime_error >()
<< "image size is different from expected: " << in.size()
<< ", " << size_;
}

OutT out{ in.size() };

constexpr auto stride_lambda = [](auto& im) -> pfft::stride_t
{
const auto [stride_x, stride_y] = im.stride();
return {static_cast<long>(stride_x), static_cast<long>(stride_y)};
};

const pfft::shape_t shape = {size_.width(), size_.height()};
const pfft::stride_t in_stride = stride_lambda(in);
const pfft::stride_t out_stride = stride_lambda(out);

pfft::c2r<ValueT>(
shape,
in_stride,
out_stride,
{ 0, 1 }, // axes
d == direction::FORWARD, // forward
reinterpret_cast<const std::complex<ValueT>*>(in.data()),
reinterpret_cast<ValueT*>(out.data()),
1.0 );

return out;
}

template < template <typename> class ImageT,
typename ContainedT,
typename ValueT = typename ContainedT::value_t,
typename OutT = image<g<ValueT>>,
typename = typename std::enable_if_t< is_imagetype_v<ImageT<ContainedT>> >
>
const cf_image& cross_correlate( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
OutT
cross_correlate( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
{
cf_image a_fft{ transform( a, direction::FORWARD ) };
cf_image b_fft{ transform( b, direction::FORWARD ) };

a_fft = b_fft * conj( a_fft );
cache().output = real( transform( a_fft, direction::REVERSE ) );
swap_quadrants( cache().output );
OutT output{ real( transform( a_fft, direction::REVERSE ) ) };
swap_quadrants( output );

return cache().output;
return output;
}

template < template <typename> class ImageT,
typename ContainedT,
typename ValueT = typename ContainedT::value_t,
typename OutT = image<g<ValueT>>,
typename = typename std::enable_if_t<
is_imagetype_v<ImageT<ContainedT>> &&
is_real_mono_pixeltype_v<ContainedT>
>
>
const cf_image& cross_correlate_real( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
OutT
cross_correlate_real( const ImageT<ContainedT>& a,
const ImageT<ContainedT>& b ) const
{
auto [a_fft, b_fft] = transform_real( a, b, direction::FORWARD );
a_fft = b_fft * conj( a_fft );
cache().output = real( transform( a_fft, direction::REVERSE ) );
swap_quadrants( cache().output );
OutT output = transform_real( a_fft, direction::REVERSE );
swap_quadrants( output );

return cache().output;
return output;
}

template < template <typename> class ImageT,
Expand Down
2 changes: 1 addition & 1 deletion openpiv/core/image.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,7 @@ class image
inline constexpr core::rect rect() const { return r_; }

// distance between pixels in x and y directions
std::tuple<size_t, size_t> stride() const
inline constexpr std::tuple<size_t, size_t> stride() const
{
return { sizeof(pixel_t), width() * sizeof(pixel_t) };
}
Expand Down
8 changes: 8 additions & 0 deletions openpiv/core/log.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,14 @@ namespace openpiv::core::logger {
{ Level::TEST, "TEST" },
} )

template <
typename Rep,
typename Period>
static std::string to_string(const std::chrono::duration<Rep, Period>& td)
{
return std::to_string(td.count());
}

template <
typename Clock,
typename Duration>
Expand Down
12 changes: 12 additions & 0 deletions openpiv/core/pixel_types.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,6 +325,18 @@ struct is_real_mono_pixeltype<g<T>> : std::true_type {};
template <typename T>
inline constexpr bool is_real_mono_pixeltype_v = is_real_mono_pixeltype<T>::value;

//
//
//
template <typename T>
struct is_complex_mono_pixeltype : std::false_type {};

template <typename T>
struct is_complex_mono_pixeltype<complex<T>> : std::true_type {};

template <typename T>
inline constexpr bool is_complex_mono_pixeltype_v = is_complex_mono_pixeltype<T>::value;


//
//
Expand Down

0 comments on commit bd58d21

Please sign in to comment.