SeparableKernel.cc

Go to the documentation of this file.

// -*- LSST-C++ -*-

/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
#include <algorithm>
#include <iterator>
#include <sstream>

#include "lsst/pex/exceptions.h"
#include "lsst/afw/math/Kernel.h"
#include "lsst/afw/math/KernelPersistenceHelper.h"
#include "lsst/afw/table/io/Persistable.cc"

namespace pexExcept = lsst::pex::exceptions;

namespace lsst {
namespace afw {

template std::shared_ptr<math::SeparableKernel> table::io::PersistableFacade<
        math::SeparableKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);

namespace math {

SeparableKernel::SeparableKernel()
        : Kernel(),
          _kernelColFunctionPtr(),
          _kernelRowFunctionPtr(),
          _localColList(0),
          _localRowList(0),
          _kernelX(0),
          _kernelY(0),
          _kernelRowCache(0),
          _kernelColCache(0) {
    _setKernelXY();
}

SeparableKernel::SeparableKernel(int width, int height, KernelFunction const& kernelColFunction,
                                 KernelFunction const& kernelRowFunction,
                                 Kernel::SpatialFunction const& spatialFunction)
        : Kernel(width, height, kernelColFunction.getNParameters() + kernelRowFunction.getNParameters(),
                 spatialFunction),
          _kernelColFunctionPtr(kernelColFunction.clone()),
          _kernelRowFunctionPtr(kernelRowFunction.clone()),
          _localColList(width),
          _localRowList(height),
          _kernelX(width),
          _kernelY(height),
          _kernelRowCache(0),
          _kernelColCache(0) {
    _setKernelXY();
}

SeparableKernel::SeparableKernel(int width, int height, KernelFunction const& kernelColFunction,
                                 KernelFunction const& kernelRowFunction,
                                 std::vector<Kernel::SpatialFunctionPtr> const& spatialFunctionList)
        : Kernel(width, height, spatialFunctionList),
          _kernelColFunctionPtr(kernelColFunction.clone()),
          _kernelRowFunctionPtr(kernelRowFunction.clone()),
          _localColList(width),
          _localRowList(height),
          _kernelX(width),
          _kernelY(height),
          _kernelRowCache(0),
          _kernelColCache(0) {
    if (kernelColFunction.getNParameters() + kernelRowFunction.getNParameters() !=
        spatialFunctionList.size()) {
        std::ostringstream os;
        os << "kernelColFunction.getNParameters() + kernelRowFunction.getNParameters() = "
           << kernelColFunction.getNParameters() << " + " << kernelRowFunction.getNParameters()
           << " != " << spatialFunctionList.size() << " = "
           << "spatialFunctionList.size()";
        throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
    }

    _setKernelXY();
}

std::shared_ptr<Kernel> SeparableKernel::clone() const {
    std::shared_ptr<Kernel> retPtr;
    if (this->isSpatiallyVarying()) {
        retPtr.reset(new SeparableKernel(this->getWidth(), this->getHeight(), *(this->_kernelColFunctionPtr),
                                         *(this->_kernelRowFunctionPtr), this->_spatialFunctionList));
    } else {
        retPtr.reset(new SeparableKernel(this->getWidth(), this->getHeight(), *(this->_kernelColFunctionPtr),
                                         *(this->_kernelRowFunctionPtr)));
    }
    retPtr->setCtr(this->getCtr());
    retPtr->computeCache(this->getCacheSize());
    return retPtr;
}

std::shared_ptr<Kernel> SeparableKernel::resized(int width, int height) const {
    std::shared_ptr<Kernel> retPtr;
    if (isSpatiallyVarying()) {
        retPtr = std::make_shared<SeparableKernel>(width, height, *_kernelColFunctionPtr,
                                                   *_kernelRowFunctionPtr, _spatialFunctionList);
    } else {
        retPtr = std::make_shared<SeparableKernel>(width, height, *_kernelColFunctionPtr,
                                                   *_kernelRowFunctionPtr);
    }
    return retPtr;
}

double SeparableKernel::computeVectors(std::vector<Pixel>& colList, std::vector<Pixel>& rowList,
                                       bool doNormalize, double x, double y) const {
    if (static_cast<int>(colList.size()) != this->getWidth() ||
        static_cast<int>(rowList.size()) != this->getHeight()) {
        std::ostringstream os;
        os << "colList.size(), rowList.size() = (" << colList.size() << ", " << rowList.size() << ") != ("
           << this->getWidth() << ", " << this->getHeight() << ") = "
           << "kernel dimensions";
        throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
    }
    if (this->isSpatiallyVarying()) {
        this->setKernelParametersFromSpatialModel(x, y);
    }

    return basicComputeVectors(colList, rowList, doNormalize);
}

SeparableKernel::KernelFunctionPtr SeparableKernel::getKernelColFunction() const {
    return _kernelColFunctionPtr->clone();
}

SeparableKernel::KernelFunctionPtr SeparableKernel::getKernelRowFunction() const {
    return _kernelRowFunctionPtr->clone();
}

std::string SeparableKernel::toString(std::string const& prefix) const {
    std::ostringstream os;
    os << prefix << "SeparableKernel:" << std::endl;
    os << prefix
       << "..x (width) function: " << (_kernelColFunctionPtr ? _kernelColFunctionPtr->toString() : "None")
       << std::endl;
    os << prefix
       << "..y (rows) function: " << (_kernelRowFunctionPtr ? _kernelRowFunctionPtr->toString() : "None")
       << std::endl;
    os << Kernel::toString(prefix + "\t");
    return os.str();
}

std::vector<double> SeparableKernel::getKernelParameters() const {
    std::vector<double> allParams = _kernelColFunctionPtr->getParameters();
    std::vector<double> yParams = _kernelRowFunctionPtr->getParameters();
    std::copy(yParams.begin(), yParams.end(), std::back_inserter(allParams));
    return allParams;
}

//
// Protected Member Functions
//

double SeparableKernel::doComputeImage(image::Image<Pixel>& image, bool doNormalize) const {
    double imSum = basicComputeVectors(_localColList, _localRowList, doNormalize);

    for (int y = 0; y != image.getHeight(); ++y) {
        image::Image<Pixel>::x_iterator imPtr = image.row_begin(y);
        for (std::vector<Pixel>::iterator colIter = _localColList.begin(); colIter != _localColList.end();
             ++colIter, ++imPtr) {
            *imPtr = (*colIter) * _localRowList[y];
        }
    }

    return imSum;
}

void SeparableKernel::setKernelParameter(unsigned int ind, double value) const {
    unsigned int const nColParams = _kernelColFunctionPtr->getNParameters();
    if (ind < nColParams) {
        _kernelColFunctionPtr->setParameter(ind, value);
    } else {
        _kernelRowFunctionPtr->setParameter(ind - nColParams, value);
    }
}

//
// Private Member Functions
//

double SeparableKernel::basicComputeVectors(std::vector<Pixel>& colList, std::vector<Pixel>& rowList,
                                            bool doNormalize) const {
    double colSum = 0.0;
    if (_kernelColCache.empty()) {
        for (unsigned int i = 0; i != colList.size(); ++i) {
            double colFuncValue = (*_kernelColFunctionPtr)(_kernelX[i]);
            colList[i] = colFuncValue;
            colSum += colFuncValue;
        }
    } else {
        int const cacheSize = _kernelColCache.size();

        int const indx = this->getKernelParameter(0) * cacheSize;

        std::vector<double>& cachedValues = _kernelColCache.at(indx);
        for (unsigned int i = 0; i != colList.size(); ++i) {
            double colFuncValue = cachedValues[i];
            colList[i] = colFuncValue;
            colSum += colFuncValue;
        }
    }

    double rowSum = 0.0;
    if (_kernelRowCache.empty()) {
        for (unsigned int i = 0; i != rowList.size(); ++i) {
            double rowFuncValue = (*_kernelRowFunctionPtr)(_kernelY[i]);
            rowList[i] = rowFuncValue;
            rowSum += rowFuncValue;
        }
    } else {
        int const cacheSize = _kernelRowCache.size();

        int const indx = this->getKernelParameter(1) * cacheSize;

        std::vector<double>& cachedValues = _kernelRowCache.at(indx);
        for (unsigned int i = 0; i != rowList.size(); ++i) {
            double rowFuncValue = cachedValues[i];
            rowList[i] = rowFuncValue;
            rowSum += rowFuncValue;

#if 0
            if (indx == cacheSize/2) {
                if (::fabs(rowFuncValue - (*_kernelRowFunctionPtr)(_kernelX[i])) > 1e-2) {
                    std::cout << indx << " " << i << " "
                              << rowFuncValue << " "
                              << (*_kernelRowFunctionPtr)(_kernelX[i])
                              << std::endl;
                }
            }
#endif
        }
    }

    double imSum = colSum * rowSum;
    if (doNormalize) {
        if ((colSum == 0) || (rowSum == 0)) {
            throw LSST_EXCEPT(pexExcept::OverflowError, "Cannot normalize; kernel sum is 0");
        }
        for (std::vector<Pixel>::iterator colIter = colList.begin(); colIter != colList.end(); ++colIter) {
            *colIter /= colSum;
        }

        for (std::vector<Pixel>::iterator rowIter = rowList.begin(); rowIter != rowList.end(); ++rowIter) {
            *rowIter /= rowSum;
        }
        imSum = 1.0;
    }
    return imSum;
}

namespace {
void _computeCache(int const cacheSize, std::vector<double> const& x,
                   SeparableKernel::KernelFunctionPtr& func, std::vector<std::vector<double> >* kernelCache) {
    if (cacheSize <= 0) {
        kernelCache->erase(kernelCache->begin(), kernelCache->end());
        return;
    }

    if (kernelCache[0].size() != x.size()) {  // invalid
        kernelCache->erase(kernelCache->begin(), kernelCache->end());
    }

    int const old_cacheSize = kernelCache->size();

    if (cacheSize == old_cacheSize) {
        return;  // nothing to do
    }

    if (cacheSize < old_cacheSize) {
        kernelCache->erase(kernelCache->begin() + cacheSize, kernelCache->end());
    } else {
        kernelCache->resize(cacheSize);
        for (int i = old_cacheSize; i != cacheSize; ++i) {
            (*kernelCache)[i].resize(x.size());
        }
    }
    //
    // Actually fill the cache
    //
    for (int i = 0; i != cacheSize; ++i) {
        func->setParameter(0, (i + 0.5) / static_cast<double>(cacheSize));
        for (unsigned int j = 0; j != x.size(); ++j) {
            (*kernelCache)[i][j] = (*func)(x[j]);
        }
    }
}
}  // namespace

void SeparableKernel::computeCache(int const cacheSize) {
    SeparableKernel::KernelFunctionPtr func;

    func = getKernelColFunction();
    _computeCache(cacheSize, _kernelY, func, &_kernelColCache);

    func = getKernelRowFunction();
    _computeCache(cacheSize, _kernelX, func, &_kernelRowCache);
}

int SeparableKernel::getCacheSize() const { return _kernelColCache.size(); };
}  // namespace math
}  // namespace afw
}  // namespace lsst