Sequential.icc

// -*- C++ -*-
// -*- coding: utf-8 -*-
//
// michael a.g. aïvázis <michael.aivazis@para-sim.com>
// parasim
// (c) 1998-2019 all rights reserved
//
 
// code guard
#if !defined(ampcor_correlators_Sequential_icc)
#error This header is an implementation detail of ampcor::correlators::Sequential
#endif
 
// interface
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
addReferenceTile(size_type pid, const constview_type & ref)
{
    // figure out the starting address of this tile in the arena
    cell_type * support = _arena + pid*(_refCells + _tgtCells);
    // adapt it into a grid
    tile_type tile(_refLayout, support);
    // move the data
    std::copy(ref.begin(), ref.end(), tile.view().begin());
    // all done
    return;
}
 
 
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
addTargetTile(size_type pid, const constview_type & tgt)
{
    // figure out the starting address of this tile in the arena
    cell_type * support = _arena + pid*(_refCells + _tgtCells) + _refCells;
    // adapt it into a grid
    tile_type tile(_tgtLayout, support);
    // move the data
    std::copy(tgt.begin(), tgt.end(), tile.view().begin());
    // all done
    return;
}
 
 
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
adjust() -> const value_type *
{
    // make a channel
    pyre::journal::debug_t channel("ampcor");
 
    // compute the dimension of the reference tiles; assume they are square
    auto refDim = _refLayout.shape()[0];
    // compute the dimension of the target tiles; assume they are square
    auto tgtDim = _tgtLayout.shape()[0];
    // compute the dimension of the correlation matrix
    auto corDim = _corLayout.shape()[0];
 
    // prelude: coarse adjustments
    auto coarseArena = _arena;
    // compute their amplitudes
    auto amplitudes = _detect(coarseArena, refDim, tgtDim);
    // adjust the reference tiles to zero mean and compute the variances
    auto refStatistics = _refStats(amplitudes, refDim, tgtDim);
    // compute the sum area tables for all possible search window placements within the target tile
    auto sat = _sat(amplitudes, refDim, tgtDim);
    // use the SATs to compute the mean amplitude of all possible window placements
    auto tgtStatistics = _tgtStats(sat, refDim, tgtDim, corDim);
    // compute the correlation hyper-surface
    auto gamma = _correlate(amplitudes, refStatistics, tgtStatistics, refDim, tgtDim, corDim);
    // find its maxima
    auto maxcor = _maxcor(gamma, corDim);
 
    // interlude: housekeeping
    delete [] gamma;
    delete [] tgtStatistics;
    delete [] sat;
    delete [] refStatistics;
    delete [] amplitudes;
 
    // refinement: refine the tiles by a factor and repeat the process with a narrower search
    // window around the location of maximum correlation
    // compute the dimension of the reference tiles; assume they are square
    auto refRefinedDim = _refRefinedLayout.shape()[0];
    // compute the dimension of the target tiles; assume they are square
    auto tgtRefinedDim = _tgtRefinedLayout.shape()[0];
    // compute the dimension of the correlation matrix
    auto corRefinedDim = _corRefinedLayout.shape()[0];
    // compute the dimension of the zoomed correlation matrix
    auto corZoomedDim = _corZoomedLayout.shape()[0];
    // ensure that the expanded target tiles where the correlation is maximum fit within the
    // search window
    _nudge(maxcor, refDim, tgtDim);
    // allocate room for the refinement area
    auto refinedArena = _refinedArena();
    // refine the reference tiles
    _refRefine(coarseArena, refinedArena);
    // collect the expanded maxcor tiles and migrate them to our new arena
    _tgtMigrate(coarseArena, maxcor, refinedArena);
    // refine the expanded target tiles in place
    _tgtRefine(refinedArena);
    // compute amplitudes
    amplitudes = _detect(refinedArena, refRefinedDim, tgtRefinedDim);
    // adjust the reference tiles to zero mean and compute the variances
    refStatistics = _refStats(amplitudes, refRefinedDim, tgtRefinedDim);
    // compute the sum area tables for all possible search window placements with the target tile
    sat = _sat(amplitudes, refRefinedDim, tgtRefinedDim);
    // use the SATs to compute the mean amplitude of all possible window placements
    tgtStatistics = _tgtStats(sat, refRefinedDim, tgtRefinedDim, corRefinedDim);
    // compute the correlation  hyper-surface
    gamma = _correlate(amplitudes, refStatistics, tgtStatistics,
                       refRefinedDim, tgtRefinedDim, corRefinedDim);
    // zoom in
    auto zoomed = _zoomcor(gamma);
    // find its maxima
    auto maxcorZoomed = _maxcor(zoomed, corZoomedDim);
    // compute the shifts and return them
    auto offsets = _offsetField(maxcorZoomed);
 
    // clean up
    delete [] maxcorZoomed;
    delete [] maxcor;
    delete [] zoomed;
    delete [] gamma;
    delete [] tgtStatistics;
    delete [] sat;
    delete [] refStatistics;
    delete [] amplitudes;
    delete [] refinedArena;
 
    // all done
    return offsets;
}
 
 
// accessors
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
arena() const -> cell_type *
{
    return _arena;
}
 
 
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
pairs() const -> size_type
{
    return _pairs;
}
 
 
// meta-methods
template <typename raster_t>
ampcor::correlators::Sequential<raster_t>::
~Sequential() {
    // release the host memory
    delete [] _arena;
    delete [] _offsets;
}
 
 
template <typename raster_t>
ampcor::correlators::Sequential<raster_t>::
Sequential(size_type pairs,
           const layout_type & refLayout, const layout_type & tgtLayout,
           size_type refineFactor, size_type refineMargin, size_type zoomFactor) :
    _pairs{ pairs },
    _refineFactor{ refineFactor },
    _refineMargin{ refineMargin },
    _zoomFactor{ zoomFactor },
    _refLayout{ refLayout },
    _tgtLayout{ tgtLayout },
    _corLayout{ tgtLayout.shape() - refLayout.shape() + index_type::fill(1) },
    _refRefinedLayout{ refineFactor * _refLayout.shape() },
    _tgtRefinedLayout{ refineFactor * (_refLayout.shape() + index_type::fill(2*refineMargin)) },
    _corRefinedLayout{ index_type::fill(2*refineFactor*refineMargin+1) },
    _corZoomedLayout { zoomFactor * _corRefinedLayout.shape() },
    _refCells{ _refLayout.size() },
    _tgtCells{ _tgtLayout.size() },
    _corCells{ _corLayout.size() },
    _refRefinedCells{ _refRefinedLayout.size() },
    _tgtRefinedCells{ _tgtRefinedLayout.size() },
    _refFootprint{ _refCells * sizeof(cell_type) },
    _tgtFootprint{ _tgtCells * sizeof(cell_type) },
    _corFootprint{ _corCells * sizeof(value_type) }, // the correlation matrix is real...
    _refRefinedFootprint{ _refRefinedCells * sizeof(cell_type) },
    _tgtRefinedFootprint{ _tgtRefinedCells * sizeof(cell_type) },
    _arena{ new cell_type[ _pairs * (_refCells+_tgtCells) ] },
    _offsets{ new value_type[ 2 * _pairs ] }
{
    // compute the footprints
    auto footprint = _pairs*(_refFootprint + _tgtFootprint);
    auto refinedFootprint = _pairs*(_refRefinedFootprint + _tgtRefinedFootprint);
    // make a channel
    pyre::journal::debug_t channel("ampcor");
    // show me
    channel
        << pyre::journal::at(__HERE__)
        << "new Sequential worker:"
        << pyre::journal::newline
        << "    pairs: " << _pairs
        << pyre::journal::newline
        << "    ref shape: " << _refLayout << ", " << _refCells << " cells"
        << pyre::journal::newline
        << "    tgt shape: " << _tgtLayout << ", " << _tgtCells << " cells"
        << pyre::journal::newline
        << "    footprint: " << (_refCells+_tgtCells) << " cells in "
        << (footprint/1024/1024) << " Mb"
        << pyre::journal::newline
        << "    refine factor: " << refineFactor
        << pyre::journal::newline
        << "    refine margin: " << refineMargin
        << pyre::journal::newline
        << "    refined ref shape: " << _refRefinedLayout << ", " << _refRefinedCells << " cells"
        << pyre::journal::newline
        << "    refined tgt shape: " << _tgtRefinedLayout << ", " << _tgtRefinedCells << " cells"
        << pyre::journal::newline
        << "    footprint: " << (_refRefinedCells+_tgtRefinedCells) << " cells in "
        << (refinedFootprint/1024/1024) << " Mb"
        << pyre::journal::newline
        << "    arena: " << _arena
        << pyre::journal::endl;
}
 
 
// debugging support
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
dump() const
{
    // dump the arena as a sequence of reference and target tiles
    pyre::journal::debug_t channel("ampcor");
 
    // sign in
    channel << pyre::journal::at(__HERE__);
    // go through all the pairs
    for (auto pid = 0; pid < _pairs; ++pid) {
        // inject the pid
        channel << "pid: " << pid << pyre::journal::newline;
        // compute the address of the reference tile in the arena
        cell_type * refLoc = _arena + pid*(_refCells + _tgtCells);
        // adapt it into a grid
        tile_type ref(_refLayout, refLoc);
        // inject it
        channel << "reference: " << pyre::journal::newline;
        for (auto idx = 0; idx < _refLayout.shape()[0]; ++idx) {
            for (auto jdx = 0; jdx < _refLayout.shape()[1]; ++jdx) {
                channel << ref[{idx, jdx}] << " ";
            }
            channel << pyre::journal::newline;
        }
 
        // compute the address of the target tile in the arena
        cell_type * tgtLoc = refLoc + _refCells;
        // adapt it into a grid
        tile_type tgt(_tgtLayout, tgtLoc);
 
        // inject it
        channel << "target: " << pyre::journal::newline;
        for (auto idx = 0; idx < _tgtLayout.shape()[0]; ++idx) {
            for (auto jdx = 0; jdx < _tgtLayout.shape()[1]; ++jdx) {
                channel << tgt[{idx, jdx}] << " ";
            }
            channel << pyre::journal::newline;
        }
 
    }
    // sing off
    channel << pyre::journal::endl;
 
    // all done
    return;
}
 
 
// implementation details: methods
 
 
 
// compute the amplitude of the signal
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_detect(const cell_type * cArena, size_type refDim, size_type tgtDim) const -> value_type *
{
    // compute the size of a reference tile
    auto refCells = refDim * refDim;
    // compute the size of a target tile
    auto tgtCells = tgtDim * tgtDim;
 
    // get a spot
    value_type * rArena = nullptr;
    // compute the number of cells whose amplitude we have to compute
    auto cells = _pairs * (refCells + tgtCells);
    // allocate room for it
    rArena = new (std::nothrow) value_type[cells]();
 
    // if something went wrong
    if (rArena == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating memory for the tile amplitudes"
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage...
    kernels::detect(cArena, cells, rArena);
 
    // all done
    return rArena;
}
 
 
 
// compute the mean and deviation of reference tiles and subtract the mean from
// each reference pixel
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_refStats(value_type * rArena, size_type refDim, size_type tgtDim) const -> value_type *
{
    // compute the size of a reference tile
    auto refCells = refDim * refDim;
    // compute the size of a target tile
    auto tgtCells = tgtDim * tgtDim;
 
    // grab a spot
    value_type * stats = nullptr;
    // allocate room for deviation: one number per reference tile
    stats = new (std::nothrow) value_type[_pairs]();
 
    // if something went wrong
    if (stats == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating memory for the variances of "
            << "the reference tiles "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage
    kernels::refStats(rArena, _pairs, refDim, refCells + tgtCells, stats);
 
    // all done
    return stats;
}
 
 
 
 
// build the sum area tables for the target tiles
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_sat(const value_type * rArena, size_type refDim, size_type tgtDim) const -> value_type *
{
 
    // compute the size of a reference tile
    auto refCells = refDim * refDim;
    // compute the size of a target tile
    auto tgtCells = tgtDim * tgtDim;
 
    // grab a spot for the sat tables
    value_type * sat = nullptr;
    // allocate memory
    sat = new (std::nothrow) value_type[_pairs*tgtCells]();
 
    // if something went wrong
    if (sat == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating memory for the sum area tables "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage
    kernels::sat(rArena, _pairs, refCells, tgtCells, tgtDim, sat);
 
    // all done
    return sat;
}
 
 
// compute the average values for all possible placements of the reference shape within the
// target tile
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_tgtStats(const value_type * dSAT,
          size_type refDim, size_type tgtDim, size_type corDim) const -> value_type *
{
    // pick a spot for the table of amplitude averages
    value_type * stats = nullptr;
    // allocate memory: one mean per placement per target tile
    stats = new (std::nothrow) value_type[_pairs*corDim*corDim];
 
    // if something went wrong
    if (stats == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating memory for the table of target amplitude averages "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage
    kernels::tgtStats(dSAT, _pairs, refDim, tgtDim, corDim, stats);
 
    // all done
    return stats;
}
 
 
// compute the correlation surface
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_correlate(const value_type * rArena,
           const value_type * refStats, const value_type * tgtSat,
           size_type refDim, size_type tgtDim, size_type corDim) const -> value_type *
{
     
    // compute the size of the reference tile
    auto refCells = refDim * refDim;
    // compute the size of the target tile
    auto tgtCells = tgtDim * tgtDim;
    // compute the size of the correlation matrix
    auto corCells = corDim * corDim;
 
    // pick a spot for the correlation matrix
    value_type * dCorrelation = nullptr;
    // compute the total number of cells in the amplitude hyper-grid
    auto size = _pairs * corCells;
    // allocate memory on the device
    dCorrelation = new (std::nothrow) value_type[size]();
 
    // if something went wrong
    if (dCorrelation == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating memory for the correlation matrix "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage
    kernels::correlate(rArena, refStats, tgtSat,
                       _pairs,
                       refCells, tgtCells, corCells, refDim, tgtDim, corDim,
                       dCorrelation);
 
    // all done
    return dCorrelation;
}
 
 
// find the locations of the correlation maxima
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_maxcor(const value_type * gamma, size_type corDim) const -> int *
{
 
    // compute the size of the correlation matrix
    auto corCells = corDim * corDim;
 
    // find a spot
    int * loc = nullptr;
    // allocate memory on the device
    loc = new (std::nothrow) int[2 * _pairs]();
    // if something went wrong
    if (loc == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "Error while allocating device memory for the location of the correlation maxima "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // engage
    kernels::maxcor(gamma, _pairs, corCells, corDim, loc);
 
    // all done
    return loc;
}
 
 
// adjust the locations of the correlation maxima so that the new target tiles fit within the
// search window
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
_nudge(int * locations, size_type refDim, size_type tgtDim) const
{
    // make sure that all locations are adjusted so that they allow enough room for the
    // {refineMargin} by moving the ULHC of the tiles so they fit
    kernels::nudge(_pairs, refDim, tgtDim, _refineMargin, locations, _offsets);
 
    // all done
    return;
}
 
 
// allocate room for the refined tiles
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_refinedArena() const -> cell_type *
{
    // grab a spot
    cell_type * arena = nullptr;
    // allocate room for it
    arena = new (std::nothrow) cell_type[_pairs * (_refRefinedFootprint + _tgtRefinedFootprint)]();
    // if something went wrong
    if (arena == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "while allocating memory for the refined tile hyper-grid "
            << pyre::journal::endl;
        // and bail
        throw std::bad_alloc();
    }
 
    // all done
    return arena;
}
 
 
// refine the reference tiles
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
_refRefine(cell_type * coarseArena, cell_type * refinedArena) const
{
    // make a channel
    pyre::journal::debug_t channel("ampcor");
 
    // get the shape the reference tile
    auto rshape = _refLayout.shape();
    // and the shape of the refined tile
    auto tshape = _refRefinedLayout.shape();
 
 
    // get number of threads. omp_get_max_threads is sometimes problematic.
    size_t nthreads=0;
    #pragma omp parallel reduction(+:nthreads)
    nthreads += 1;
 
    // step 1: initiate FFT processor 
    isce3::signal::Signal<float> * procFFT = nullptr;
    procFFT = new (std::nothrow) isce3::signal::Signal<float>(nthreads);
 
    // if something went wrong
    if (procFFT == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "while instanciating a isce3::signal::Signal FFT processor" 
            << pyre::journal::endl;
        // and bail
        throw std::runtime_error("while instanciating a isce3::signal::Signal FFT processor");
    }
 
 
    // step 2: forward FFT from {coarseArena} to {refinedArena}
    // the plan characteristics
    int dim = 2;
    int fwdRanks[] = { static_cast<int>(rshape[0]), static_cast<int>(rshape[1]) };
    //  the data layout of the coarse tiles
    // the tile shape stays the same
    int fwdIEmbed[] = { static_cast<int>(rshape[0]), static_cast<int>(rshape[1]) };
    // the data is densely packed
    int fwdIStride = 1;
    // the distance between reference tiles
    int fwdIDist = _refCells + _tgtCells;
    // the data layout of the refined tiles
    // the tile shape stays the same
    int fwdOEmbed[] = { static_cast<int>(tshape[0]), static_cast<int>(tshape[1]) };
    // the data is densely packed
    int fwdOStride = 1;
    // the distance between reference tiles in the refined arena
    int fwdODist = _refRefinedCells + _tgtRefinedCells;
 
    // Set the FFT forward plan for upsampling reference tiles
    // This involves out-of-place FFT 
    procFFT->fftPlanForward(reinterpret_cast<std::complex<float> *>(coarseArena),
                            reinterpret_cast<std::complex<float> *>(refinedArena),
                            dim, fwdRanks, _pairs,
                            fwdIEmbed, fwdIStride, fwdIDist,
                            fwdOEmbed, fwdOStride, fwdODist, -1);
 
 
 
 
    // step 3: inverse FFT of the tiles back into the refined arena
    // the plan characteristics
    int revRanks[] = { static_cast<int>(tshape[0]), static_cast<int>(tshape[1]) };
    // the data layout of the transformed reference tiles
    int revIEmbed[] = { static_cast<int>(tshape[0]), static_cast<int>(tshape[1]) };
    // the data is densely packed
    int revIStride = 1;
    // the distance between reference tiles
    int revIDist = _refRefinedCells + _tgtRefinedCells;
    // the inverse FFT tiles have identical layout
    int revOEmbed[] = { static_cast<int>(tshape[0]), static_cast<int>(tshape[1]) };
    // the data is densely packed
    int revOStride = 1;
    // the distance between reference tiles
    int revODist = _refRefinedCells + _tgtRefinedCells;
 
    procFFT->fftPlanBackward(reinterpret_cast<std::complex<float> *>(refinedArena),
                             reinterpret_cast<std::complex<float> *>(refinedArena),
                             //dummy,
                             dim, revRanks, _pairs,
                             revIEmbed, revIStride, revIDist,
                             revOEmbed, revOStride, revODist, 1);
 
 
    // step 4: Actual refinement (upsampling) of the reference tiles:
    // FFT -> frequencies shuffling -> back FFT
    procFFT->upsample2D(reinterpret_cast<std::complex<float> *>(coarseArena),
                        reinterpret_cast<std::complex<float> *>(refinedArena), 
                        _refineFactor);
 
    // cleaning
    delete procFFT;
 
    // all done
    return;
 
}
 
 
 
 
 
// migrate the expanded unrefined target tiles into the {refinedArena}
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
_tgtMigrate(cell_type * coarseArena, int * locations, cell_type * refinedArena) const
{
    // make a channel
    pyre::journal::debug_t channel("ampcor");
 
    // the reference tile shape
    auto refShape = _refLayout.shape();
    // the target tile shape
    auto tgtShape = _tgtLayout.shape();
    // the refined reference tile shape
    auto refRefinedShape = _refRefinedLayout.shape();
    // the refined target tile shape
    auto tgtRefinedShape = _tgtRefinedLayout.shape();
 
    // unpack the dimensions
    auto refDim = refShape[0];
    auto tgtDim = tgtShape[0];
    auto refRefinedDim = refRefinedShape[0];
    auto tgtRefinedDim = tgtRefinedShape[0];
 
    // compute the dimension of the expanded maxcor tile
    auto expDim = refDim + 2 * _refineMargin;
 
    // engage...
    kernels::migrate(coarseArena, _pairs,
                     refDim, tgtDim, expDim,
                     refRefinedDim, tgtRefinedDim,
                     locations,
                     refinedArena);
 
    // all done
    return;
}
 
 
// refine the target tiles around the locations of the correlation maxima
template <typename raster_t>
void
ampcor::correlators::Sequential<raster_t>::
_tgtRefine(cell_type * refinedArena) const
{
    // N.B.: assumes {locations} are already nudged and on the CPU...
    // make a channel
    pyre::journal::debug_t channel("ampcor");
 
    // the shape of the target tile
    //auto tgtShape = _tgtLayout.shape();
    // the shape the refined target tile
    auto tgtRefShape = _tgtRefinedLayout.shape();
    // the shape of the expanded target tile
    auto expShape = _refLayout.shape() + index_type::fill(2*_refineMargin);
 
    // N.B.: the expanded maxcor target tiles are expected to have already been moved to the
    // refined arena after the maxcor locations were determined
 
 
    // get number of threads. omp_get_max_threads is sometimes problematic.
    size_t nthreads=0;
    #pragma omp parallel reduction(+:nthreads)
    nthreads += 1;
 
    // step 1: initiate FFT processor 
    isce3::signal::Signal<float> * procFFT = nullptr;
    procFFT = new (std::nothrow) isce3::signal::Signal<float>(nthreads);
 
    // if something went wrong
    if (procFFT == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "while instanciating a isce3::signal::Signal FFT processor" 
            << pyre::journal::endl;
        // and bail
        throw std::runtime_error("while instanciating a isce3::signal::Signal FFT processor");
    }
 
    
 
    // step 2: forward FFT in place in {refinedArena}
    // the plan characteristics
    int dim = 2;
    // use the shape of the expanded target tile
    int fwdRanks[] = { static_cast<int>(expShape[0]), static_cast<int>(expShape[1]) };
    // this tile is already occupying its destination location in {refinedArena}
    int fwdIEmbed[] = { static_cast<int>(tgtRefShape[0]), static_cast<int>(tgtRefShape[1]) };
    // the data is dense
    int fwdIStride = 1;
    // the distance between tiles
    int fwdIDist = _refRefinedCells + _tgtRefinedCells;
    // the destination of the forward FFT has indentical layout
    int fwdOEmbed[] = { static_cast<int>(tgtRefShape[0]), static_cast<int>(tgtRefShape[1]) };
    // the data is dense
    int fwdOStride = 1;
    // the distance between tiles
    int fwdODist = _refRefinedCells + _tgtRefinedCells;
 
    // the address of the first expanded target tile
    auto firstTile = reinterpret_cast<std::complex<float> *>(refinedArena + _refRefinedCells);
 
    procFFT->fftPlanForward(firstTile, firstTile,
                            dim, fwdRanks, _pairs,
                            fwdIEmbed, fwdIStride, fwdIDist,
                            fwdOEmbed, fwdOStride, fwdODist, -1);
 
 
 
    // step 3: inverse FFT of the refined tiles back into the refined arena
    // the plan characteristics
    int revRanks[] = { static_cast<int>(tgtRefShape[0]), static_cast<int>(tgtRefShape[1]) };
    // the data layout of the transformed reference tiles
    int revIEmbed[] = { static_cast<int>(tgtRefShape[0]), static_cast<int>(tgtRefShape[1]) };
    // the data is densely packed
    int revIStride = 1;
    // the distance between reference tiles
    int revIDist = _refRefinedCells + _tgtRefinedCells;
    // the inverse FFT tiles have identical layout
    int revOEmbed[] = { static_cast<int>(tgtRefShape[0]), static_cast<int>(tgtRefShape[1]) };
    // the data is densely packed
    int revOStride = 1;
    // the distance between reference tiles
    int revODist = _refRefinedCells + _tgtRefinedCells;
 
    procFFT->fftPlanBackward(firstTile, firstTile,
                             dim, revRanks, _pairs,
                             revIEmbed, revIStride, revIDist,
                             revOEmbed, revOStride, revODist, 1);
 
 
    // step4: Upsample!
    procFFT->upsample2D(firstTile, firstTile, _refineFactor);
 
    
    // no need of plan anymore
    delete procFFT;
 
    // all done
    return;
 
}
 
 
 
 
 
 
 
 
// zoom the correlation matrix
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_zoomcor(value_type * gamma) const -> value_type *
{
    // make a channel
    pyre::journal::debug_t channel("ampcor");
 
    // get the shape of the incoming correlation matrix
    auto corShape = _corRefinedLayout.shape();
    // extract the dimension
    int corDim = corShape[0];
    // get the shape of the zoomed correlation matrix
    auto zmdShape = _corZoomedLayout.shape();
    // compute the number if cells in each zoomed correlation matrix
    auto zmdCells = _corZoomedLayout.size();
    // extract the dimension
    int zmdDim = zmdShape[0];
 
    // step 1: up-cast and embed
    auto scratch = kernels::r2c(gamma, _pairs, corDim, zmdDim);
 
 
    // step 2: forward FFT from the incoming gamma to the zoomed gamma
    // the plan characteristics
    int dim = 2;
    int ranks[] = { corDim, corDim };
    // the input layout
    int iEmbed[] = { zmdDim, zmdDim };
    // the data is densely packed
    int iStride = 1;
    // the distance between correlation matrices
    int iDist = zmdCells;
    // the output layout
    int oEmbed[] = { zmdDim, zmdDim };
    // the data is densely packed
    int oStride = 1;
    // the distance between successive correlation matrices
    int oDist = zmdCells;
 
    // get number of threads. omp_get_max_threads is sometimes problematic.
    size_t nthreads=0;
    #pragma omp parallel reduction(+:nthreads)
    nthreads += 1;
 
    // step 3: initiate FFT processor 
    isce3::signal::Signal<float> * procFFT = nullptr;
    procFFT = new (std::nothrow) isce3::signal::Signal<float>(nthreads);
 
    // if something went wrong
    if (procFFT == nullptr) {
        // make a channel
        pyre::journal::error_t error("ampcor");
        // complain
        error
            << pyre::journal::at(__HERE__)
            << "while instanciating a isce3::signal::Signal FFT processor" 
            << pyre::journal::endl;
        // and bail
        throw std::runtime_error("while instanciating a isce3::signal::Signal FFT processor");
    }
 
 
    // Forward FFT plan - in-place transformation
    procFFT->fftPlanForward(scratch, scratch,
                            dim, ranks, _pairs,
                            iEmbed, iStride, iDist,
                            oEmbed, oStride, oDist, -1);
 
    // Reverse FFT plan - in-place transformation
    ranks[0] = zmdDim;
    ranks[1] = zmdDim;
    procFFT->fftPlanBackward(scratch, scratch,
                             dim, ranks, _pairs,
                             oEmbed, oStride, oDist,
                             oEmbed, oStride, oDist, 1);
 
    // Upsample the correlation map
    // fwd FFT then spectrum shuffling then rev FFT
    procFFT->upsample2D(scratch, scratch, _zoomFactor);
 
 
    // no need of plan anymore
    delete procFFT;
 
    // Convert complex to real
    auto zoomed = kernels::c2r(scratch, _pairs, zmdDim);
 
    // clean up
    delete scratch;
 
    // all done
    return zoomed;
}
 
 
 
 
// assemble the offsets
template <typename raster_t>
auto
ampcor::correlators::Sequential<raster_t>::
_offsetField(const int * fine) -> const value_type *
{
    // compute the search margin
    auto margin = (_tgtLayout.shape()[0] - _refLayout.shape()[0]) / 2;
    // compute the overall zoom factor
    auto zoom = _refineFactor * _zoomFactor;
 
    // launch the kernel that does the work
    kernels::offsetField(fine, _pairs, margin, _refineMargin, zoom, _offsets);
 
    // all done
    return _offsets;
}
 
 
 
// end of file