-
Christof Kaufmann authoredChristof Kaufmann authored
image.cpp 100.26 KiB
#include <algorithm>
#include <numeric>
#include <type_traits>
#include <sstream>
#include "filesystem.h"
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <gdal_utils.h>
#include <cpl_string.h>
#include "logging.h"
#include "image.h"
#include "geoinfo.h"
namespace {
template <typename Imgtype>
imagefusion::Image merge_impl(std::vector<Imgtype> const& images) {
std::vector<cv::Mat> cvImages;
cvImages.reserve(images.size());
for (auto const& i : images)
cvImages.push_back(i.cvMat()); // OpenCV shared copy
// unsigned int totalChannels = std::accumulate(images.begin(), images.end(), 0, [] (Imgtype const& i, unsigned int count) { return count + i.channels(); });
imagefusion::Image multichannel;
cv::merge(cvImages, multichannel.cvMat());
return multichannel;
}
std::string getGDALMemOptionString(imagefusion::ConstImage const& i) {
char szPtrValue[128] = { '\0' };
int nRet = CPLPrintPointer(szPtrValue,
reinterpret_cast<void*>(const_cast<uchar*>(i.cvMat().ptr())),
sizeof(szPtrValue));
szPtrValue[nRet] = 0;
std::ostringstream ss;
ss << "MEM:::DATAPOINTER=" << szPtrValue
<< ",PIXELS=" << i.width()
<< ",LINES=" << i.height()
<< ",BANDS=" << i.channels()
<< ",DATATYPE=" << GDALGetDataTypeName(toGDALDepth(i.type()))
<< ",PIXELOFFSET=" << i.cvMat().elemSize()
<< ",LINEOFFSET=" << i.cvMat().ptr(1) - i.cvMat().ptr(0)
<< ",BANDOFFSET=" << i.cvMat().elemSize1();
return ss.str();
}
} /* anonymous namespace */
namespace imagefusion {
// Implementation
GDALDataset const* ConstImage::asGDALDataset() const {
// This reads the OpenCV Mat as GDAL Dataset. See ConstImage::write for more info!
std::string optStr = getGDALMemOptionString(*this);
return (GDALDataset*) GDALOpen(optStr.c_str(), GA_ReadOnly);
}
GDALDataset* Image::asGDALDataset() {
// This reads the OpenCV Mat as GDAL Dataset. See ConstImage::write for more info!
std::string optStr = getGDALMemOptionString(*this);
return (GDALDataset*) GDALOpen(optStr.c_str(), GA_Update);
}
ConstImage const& ConstImage::write(std::string const& filename, GeoInfo const& gi, FileFormat format) const {
GDALAllRegister();
if (format == FileFormat::unsupported) {
fs::path p = filename;
std::string ext = p.extension().string();
format = FileFormat::fromFileExtension(ext);
if (format == FileFormat::unsupported) {
IF_THROW_EXCEPTION(file_format_error("Cannot auto detect image format for file extension " + ext +
". Please specify image format explicitly or read it from an input image file!"))
<< boost::errinfo_file_name(filename);
}
}
try {
if (format == FileFormat("GTiff"))
return write(filename, to_string(format), {std::make_pair<std::string,std::string>("COMPRESS", "LZW")}, gi);
else
return write(filename, to_string(format), {}, gi);
}
catch(boost::exception& ex) {
ex << errinfo_file_format(to_string(format));
throw;
}
}
ConstImage const& ConstImage::write(std::string const& filename, std::string const& driverName, std::vector<std::pair<std::string,std::string>> const& options_vec, GeoInfo const& gi) const {
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName(driverName.c_str());
if (!driver) {
IF_THROW_EXCEPTION(file_format_error("GDAL driver for image file format " + driverName + " not available! "
"Maybe you have to recompile GDAL with more dependencies."))
<< boost::errinfo_file_name(filename);
}
// To write as png or jpg, we need to use CreateCopy method, from a memory ("MEM") Dataset:
// https://gis.stackexchange.com/questions/134000/gdal-api-cant-save-image-in-some-formats
// To avoid making an additional copy, we use the memory of the OpenCV Mat for the GDAL Dataset.
// For that we need to convert some options (like pointer to image memory) to a string and use it as filename:
// https://gis.stackexchange.com/questions/196048/how-to-reuse-memory-pointer-of-gdal-memory-driver
// https://lists.osgeo.org/pipermail/gdal-dev/2016-June/044556.html
GDALDataset* poSrcDS = const_cast<GDALDataset*>(asGDALDataset());
// add GeoInfo first time (important for color table)
gi.addTo(poSrcDS);
// write to file
char **options = nullptr;
for (auto const& p : options_vec)
options = CSLSetNameValue(options, p.first.c_str(), p.second.c_str());
GDALDataset* poDstDS = driver->CreateCopy(filename.c_str(), poSrcDS, FALSE, options, nullptr, nullptr);
// add GeoInfo second time (important for custom metadata domains)
gi.addTo(poDstDS);
bool shouldHaveWrittenColorTable = !gi.colorTable.empty() && type() == Type::uint8x1;
GDALClose(poSrcDS);
if (options)
CSLDestroy(options);
if (!poDstDS) {
IF_THROW_EXCEPTION(runtime_error("Could not create the output file " + filename + " with driver " + driverName +
" (" + FileFormat(driverName).longName() + ") and type " + to_string(type()) + "."))
<< errinfo_file_format(driverName)
<< errinfo_image_type(type())
<< boost::errinfo_file_name(filename);
}
GDALClose(poDstDS);
// check for changed color table and warn if (shouldHaveWrittenColorTable) {
GeoInfo test{filename};
gi.compareColorTables(test, false);
}
return *this;
}
Image& Image::read(std::string const& filename,
std::vector<int> channels /* {} means all channels */,
Rectangle r /* {0, 0, 0, 0} means whole image */,
bool flipH, bool flipV, bool ignoreColorTable, InterpMethod interp)
{
GDALAllRegister();
// open file
FileFormat format = FileFormat::fromFile(filename);
if (format == FileFormat::unsupported)
IF_THROW_EXCEPTION(runtime_error("Could not open image '" + filename + "' with GDAL to read the image. "
"Either the file does not exist or GDAL could not find an appropriate driver to read the image."))
<< boost::errinfo_file_name(filename);
// open dataset
GeoInfo gi{filename};
unsigned int rastercount = gi.channels;
GDALDataset* gdal_img = nullptr;
if (gi.hasSubdatasets()) {
gdal_img = gi.openVrtGdalDataset(channels, interp);
rastercount = gdal_img->GetRasterCount();
GeoInfo gi_sub;
gi_sub.readFrom(gdal_img);
if (gi_sub.baseType == Type::invalid) {
std::string types;
for (unsigned int i = 1; i <= rastercount; ++i) {
GDALDataType gdal_type = gdal_img->GetRasterBand(i)->GetRasterDataType();
types += to_string(toBaseType(gdal_type)) + ", ";
}
types.pop_back(); types.pop_back();
GDALClose(gdal_img);
IF_THROW_EXCEPTION(image_type_error("The selected subdatasets have different data types: " + types + ". Cannot import image with different types."));
}
channels.clear(); // to get all channels in virtual dataset
}
else {
// check acquired channels
if (!channels.empty()) {
for (int& c : channels) { // in GDAL channels are called bands and are 1 based
if (c >= static_cast<int>(rastercount) || c < 0)
IF_THROW_EXCEPTION(image_type_error("You acquired a channel (" + std::to_string(c)
+ ") that does not exist. The image only has "
+ std::to_string(rastercount) + " channels."))
<< boost::errinfo_file_name(filename);
++c;
}
rastercount = channels.size();
}
gdal_img = (GDALDataset*) GDALOpen(filename.c_str(), GA_ReadOnly);
}
if (!gdal_img)
IF_THROW_EXCEPTION(file_format_error("Could not open image " + filename + " for unknown reasons. This might be a bug."))
<< boost::errinfo_file_name(filename);
// check acquired region
int cols = gdal_img->GetRasterXSize();
int rows = gdal_img->GetRasterYSize();
if (r.width == 0)
r.width = cols - r.x;
if (r.height == 0)
r.height = rows - r.y;
if (r.x + r.width > cols || r.y + r.height > rows // check bounds
|| r.x < 0 || r.y < 0 || r.width < 0 || r.height < 0) // check negative values
{
IF_THROW_EXCEPTION(size_error("The region you acquired " + to_string(r) + " has negative size or "
"goes out of bounds of the image with size ("
+ std::to_string(cols) + " x " + std::to_string(rows) + ")."))
<< boost::errinfo_file_name(filename);
}
// create image memory
GDALDataType gdal_type = gdal_img->GetRasterBand(1)->GetRasterDataType();
int cv_type = CV_MAKETYPE(toCVType(toBaseType(gdal_type)), rastercount);
img.create(r.height, r.width, cv_type);
// read
uint8_t* pData = img.ptr();
if (flipH)
pData += (r.width - 1) * img.elemSize();
if (flipV)
pData += (r.height - 1) * (img.ptr(1) - img.ptr(0));
auto ret = gdal_img->RasterIO(
GF_Read,
r.x, /* nXOff */
r.y, /* nYOff */
r.width, /* nXSize */
r.height, /* nYSize */
pData, /* pData */
r.width, /* nBufXSize */
r.height, /* nBufYSize */
gdal_type, /* eBufType */
rastercount, /* nBandCount */
channels.empty() ? nullptr : channels.data(), /* panBandMap, nullptr means all bands */
img.elemSize() * (flipH ? -1 : 1), /* nPixelSpace, how many bytes to the next value of the same channel*/
(img.ptr(1) - img.ptr(0)) * (flipV ? -1 : 1), /* nLineSpace */
img.elemSize1() /* nBandSpace */
// NULL /* psExtraArg */
);
if (ret == CE_Failure) {
GDALClose(gdal_img);
IF_THROW_EXCEPTION(runtime_error("Could not read from image " + filename +
" of type " + to_string(type()) + "."))
<< errinfo_image_type(type())
<< boost::errinfo_file_name(filename);
}
// check for color table (indexed colors) and replace image with converted colors image
GDALColorTable* ct = gdal_img->GetRasterBand(1)->GetColorTable();
if (ct && !ignoreColorTable) {
if (ct->GetPaletteInterpretation() == GPI_CMYK || ct->GetPaletteInterpretation() == GPI_HLS || (ct->GetPaletteInterpretation() == GPI_RGB && rastercount > 1))
IF_THROW_EXCEPTION(file_format_error("Color indexed images with RGB interpretation are only supported for single-index images. CMYK and HLS are not supported at all."))
<< boost::errinfo_file_name(filename);
if (gdal_type != GDT_Byte)
IF_THROW_EXCEPTION(file_format_error("This is a color indexed image, which is not of type uint8 (or GDT_Byte). I did not know this is possible. Please help me to fix this!"))
<< boost::errinfo_file_name(filename);
bool isGray = true;
bool hasAlpha = false;
int ctsize = ct->GetColorEntryCount();
if (ct->GetPaletteInterpretation() == GPI_RGB) {
// test only the occuring colors to check for the number of required target channels
for (int y = 0; y < r.height && (isGray || !hasAlpha); ++y) { for (int x = 0; x < r.width && (isGray || !hasAlpha); ++x) {
int i = at<uint8_t>(x, y, 0);
if (i >= ctsize)
continue;
GDALColorEntry const* ce = ct->GetColorEntry(i);
if (ce->c1 > 255 || ce->c2 > 255 || ce->c3 > 255 || ce->c4 > 255)
IF_THROW_EXCEPTION(file_format_error("The color index has values that do not fit into 8 bit. I did not know this is possible. Please help me to fix this!"))
<< boost::errinfo_file_name(filename);
if (ce->c1 != ce->c2 || ce->c2 != ce->c3)
isGray = false;
if (ce->c4 != 255)
hasAlpha = true;
}
}
}
int targetChannels = isGray ? 1 : 3;
if (hasAlpha)
++targetChannels;
if (ct->GetPaletteInterpretation() == GPI_RGB) {
SPDLOG_SINGLE_INFO("Note, interpreted the RGBA color table of {} as {}{} This results in {} channel(s) in the converted image.",
filename,
(isGray ? "Gray scale (since all used entries are gray) " : "RGB (since there are non-gray colors) "),
(hasAlpha ? "with Alpha channel (since the alpha values vary)." : "without Alpha channel (since the alpha values are all 255)."),
targetChannels);
}
cv_type = CV_MAKETYPE(toCVType(Type::uint8), targetChannels);
cv::Mat indexed = img;
// std::cout << "indexed:\n" << img << std::endl;
img.create(r.height, r.width, cv_type);
GDALColorEntry black{0, 0, 0, 255};
for (int y = 0; y < r.height; ++y) {
for (int x = 0; x < r.width; ++x) {
int i = indexed.at<uint8_t>(/*cv order*/ y, x);
GDALColorEntry const* ce;
if (i >= ctsize) {
std::cout << "WARNING: Cannot convert the indexed values at (" << x << ", " << y << "), because it is too large: "
<< i << " >= " << ctsize << " = table size. (Maybe you read out of bounds?) Using black instead." << std::endl;
ce = &black;
}
else
ce = ct->GetColorEntry(i);
switch(targetChannels) {
case 4:
at<uint8_t>(x, y, 3) = cv::saturate_cast<uint8_t>(ce->c4);
[[fallthrough]];
case 3:
at<uint8_t>(x, y, 2) = cv::saturate_cast<uint8_t>(ce->c3);
at<uint8_t>(x, y, 1) = cv::saturate_cast<uint8_t>(ce->c2);
at<uint8_t>(x, y, 0) = cv::saturate_cast<uint8_t>(ce->c1);
break;
case 2:
at<uint8_t>(x, y, 1) = cv::saturate_cast<uint8_t>(ce->c4); // gray-alpha
[[fallthrough]];
case 1:
at<uint8_t>(x, y, 0) = cv::saturate_cast<uint8_t>(ce->c1);
}
}
}
// std::cout << "converted:\n" << img << std::endl;
}
GDALClose(gdal_img);
return *this;
}
Rectangle ConstImage::getCropWindow() const {
cv::Size wholeSize;
cv::Point offset;
img.locateROI(wholeSize, offset);
return Rectangle{offset.x, offset.y, width(), height()};
}
Size ConstImage::getOriginalSize() const {
cv::Size wholeSize;
cv::Point offset;
img.locateROI(wholeSize, offset);
return Size{wholeSize.width, wholeSize.height};
}
ConstImage& ConstImage::crop(Rectangle r) {
ConstImage temp = sharedCopy(r);
img = std::move(temp.img);
return *this;
}
ConstImage& ConstImage::uncrop() {
Size s = getOriginalSize();
return adjustCropBorders(s.height, s.height, s.width, s.width);
}
ConstImage& ConstImage::adjustCropBorders(int extendTop, int extendBottom, int extendLeft, int extendRight) {
// check for resulting zero size
cv::Size wholeSize;
cv::Point offset;
img.locateROI(wholeSize, offset);
int newLeft = offset.x - extendLeft;
int newRight = offset.x + width() + extendRight;
int newTop = offset.y - extendTop;
int newBottom = offset.y + height() + extendBottom;
newLeft = std::max(newLeft, 0);
newRight = std::min(newRight, wholeSize.width);
newTop = std::max(newTop, 0);
newBottom = std::min(newBottom, wholeSize.height);
const int newWidth = newRight - newLeft;
const int newHeight = newBottom - newTop;
if (newWidth <= 0 || newHeight <= 0)
IF_THROW_EXCEPTION(size_error("Adjusting crop borders would result in zero size (here, width: " + std::to_string(newWidth)
+ ", height: " + std::to_string(newHeight) + "), which is not supported."));
// non-zero size, do it!
img.adjustROI(extendTop, extendBottom, extendLeft, extendRight);
return *this;
}
Image& Image::copyValuesFrom(ConstImage const& src, ConstImage const& mask) {
assert(src.width() == width());
assert(src.height() == height());
assert(src.channels() == channels());
assert(src.type() == type());
if (mask.empty())
src.img.copyTo(img);
else
src.img.copyTo(img, mask.img);
return *this;
}
Image ConstImage::clone() const { if (empty())
return {};
cv::Size wholeSize;
cv::Point offset;
img.locateROI(wholeSize, offset);
Image ret;
if (wholeSize.height == height() && wholeSize.width == width()) {
// not cropped, just clone
ret.img = img.clone();
return ret;
}
// this is cropped, so uncrop it first and then clone and recrop
Image temp;
temp.img = img;
temp.uncrop();
ret.img = temp.img.clone();
ret.crop(Rectangle{offset.x, offset.y, width(), height()});
return ret;
}
Image ConstImage::clone(Rectangle const& r) const {
ConstImage temp = sharedCopy(r);
Image ret;
ret.img = temp.img.clone();
return ret;
}
namespace {
struct BilinearFilterHelper {
ConstImage const& img_tl;
ConstImage const& img_tr;
ConstImage const& img_bl;
ConstImage const& img_br;
double dx;
double dy;
Size s;
Type fullType;
template<Type baseType>
Image operator()() const {
using type = typename DataType<baseType>::base_type;
int chans = getChannels(fullType);
Image cropped{s, fullType};
double rounding = isIntegerType(baseType) ? 0.5 : 0;
double weight_tl = (1 - dx) * (1 - dy);
double weight_tr = dx * (1 - dy);
double weight_bl = (1 - dx) * dy ;
double weight_br = dx * dy ;
if (dx == 0) {
// only interpolate y
for (int y = 0; y < s.height; ++y) {
for (int x = 0; x < s.width; ++x) {
for (int c = 0; c < chans; ++c) {
double tl = img_tl.at<type>(x, y, c) * weight_tl;
double bl = img_bl.at<type>(x, y, c) * weight_bl;
cropped.at<type>(x, y, c) = static_cast<type>(tl + bl + rounding);
}
}
}
}
else if (dy == 0) {
// only interpolate x
for (int y = 0; y < s.height; ++y) { for (int x = 0; x < s.width; ++x) {
for (int c = 0; c < chans; ++c) {
double tl = img_tl.at<type>(x, y, c) * weight_tl;
double tr = img_tr.at<type>(x, y, c) * weight_tr;
cropped.at<type>(x, y, c) = static_cast<type>(tl + tr + rounding);
}
}
}
}
else {
// interpolate x and y
for (int y = 0; y < s.height; ++y) {
for (int x = 0; x < s.width; ++x) {
for (int c = 0; c < chans; ++c) {
double tl = img_tl.at<type>(x, y, c) * weight_tl;
double tr = img_tr.at<type>(x, y, c) * weight_tr;
double bl = img_bl.at<type>(x, y, c) * weight_bl;
double br = img_br.at<type>(x, y, c) * weight_br;
cropped.at<type>(x, y, c) = static_cast<type>(tl + tr + bl + br + rounding);
}
}
}
}
return cropped;
}
};
}
Image ConstImage::clone(Coordinate const& topleft, Size s) const {
double dx = topleft.x - std::floor(topleft.x);
double dy = topleft.y - std::floor(topleft.y);
if (dx == 0 && dy == 0)
return clone(Rectangle{static_cast<int>(topleft.x), static_cast<int>(topleft.y), s.width, s.height});
int x0 = static_cast<int>(std::floor(topleft.x));
int x1 = static_cast<int>(std::ceil(topleft.x));
int y0 = static_cast<int>(std::floor(topleft.y));
int y1 = static_cast<int>(std::ceil(topleft.y));
ConstImage img_tl = sharedCopy(Rectangle{x0, y0, s.width, s.height});
ConstImage img_tr = sharedCopy(Rectangle{x1, y0, s.width, s.height});
ConstImage img_bl = sharedCopy(Rectangle{x0, y1, s.width, s.height});
ConstImage img_br = sharedCopy(Rectangle{x1, y1, s.width, s.height});
return CallBaseTypeFunctor::run(BilinearFilterHelper{img_tl, img_tr, img_bl, img_br, dx, dy, s, type()}, basetype());
// this code uses much more memory than the above functor, since it generates a temporary double image and image-wide operations. It is only slightly faster.
// cv::Mat temp;
// cv::addWeighted(img_tl.cvMat(), weight_tl, img_tr.cvMat(), weight_tr, 0, temp, CV_64F);
// cv::addWeighted(temp, 1, img_bl.cvMat(), weight_bl, 0, temp, CV_64F);
// cv::addWeighted(temp, 1, img_br.cvMat(), weight_br, 0, temp, CV_64F);
// Image cropped;
// temp.convertTo(cropped.cvMat(), toCVType(basetype())); // converting from double to integer type will also round the result
// return cropped;
}
Image ConstImage::warp(GeoInfo const& from, GeoInfo const& to, InterpMethod method) const {
if (!from.hasGeotransform() || !to.hasGeotransform())
IF_THROW_EXCEPTION(invalid_argument_error("For warping you need to specify from which projection system to which should be warped. "
"One of your arguments does not specify any geotransformation (GCPs are currently ignored)."));
GDALDataset* poSrcDS = const_cast<GDALDataset*>(asGDALDataset());
from.addTo(poSrcDS);
Size sz = to.size;
if (sz.width <= 0 || sz.height <= 0) {
// estimate output size from source
char* projStr;
to.geotransSRS.exportToWkt(&projStr);
if (std::strlen(projStr) <= 0)
IF_THROW_EXCEPTION(logic_error("no reference system"));
GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;
void* pTransformArg = GDALCreateGenImgProjTransformer(poSrcDS,
poSrcDS->GetProjectionRef(),
nullptr,
projStr,
FALSE, 0.0, 1);
double geoTransformOut[6];
int widthOut;
int heightOut;
double extentsOut[4]; // xmin, ymin, xmax, ymax in projection space
#ifdef DEBUG
CPLErr eErr =
#endif
GDALSuggestedWarpOutput2(poSrcDS,
pfnTransformer,
pTransformArg,
geoTransformOut,
&widthOut,
&heightOut,
extentsOut,
0);
GDALDestroyGenImgProjTransformer(pTransformArg);
CPLFree(projStr);
CPLAssert(eErr == CE_None);
// std::cout << "suggested... width: " << widthOut << ", height: " << heightOut
// << ", extents: (" << extentsOut[0] << ", " << extentsOut[1] << "), (" << extentsOut[2] << ", " << extentsOut[3] << ")" << std::endl;
Coordinate c1 = to.geotrans.projToImg({extentsOut[0], extentsOut[1]});
Coordinate c2 = to.geotrans.projToImg({extentsOut[2], extentsOut[3]});
sz.width = static_cast<int>(std::ceil(std::max(c1.x, c2.x)));
sz.height = static_cast<int>(std::ceil(std::max(c1.y, c2.y)));
}
Image warped{sz, type()};
GDALDataset* poDstDS = warped.asGDALDataset();
to.addTo(poDstDS);
// Setup warp options.
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
if (method == InterpMethod::nearest)
psWarpOptions->eResampleAlg = GRA_NearestNeighbour; // for data images
else if (method == InterpMethod::bilinear)
psWarpOptions->eResampleAlg = GRA_Bilinear;
else if (method == InterpMethod::cubic)
psWarpOptions->eResampleAlg = GRA_Cubic;
else // method == InterpMethod::cubicspline
psWarpOptions->eResampleAlg = GRA_CubicSpline;
psWarpOptions->hSrcDS = poSrcDS;
psWarpOptions->hDstDS = poDstDS;
psWarpOptions->nBandCount = channels(); // required for nodata values
psWarpOptions->panSrcBands = (int*) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
psWarpOptions->panDstBands = (int*) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
for (unsigned int i = 0; i < channels(); ++i) {
psWarpOptions->panSrcBands[i] = i + 1;
psWarpOptions->panDstBands[i] = i + 1;
}
// nodata values
if (from.hasNodataValue()) {
psWarpOptions->padfSrcNoDataReal = (double*) CPLMalloc(channels() * sizeof(double));
psWarpOptions->padfSrcNoDataImag = (double*) CPLMalloc(channels() * sizeof(double)); psWarpOptions->padfDstNoDataReal = (double*) CPLMalloc(channels() * sizeof(double));
psWarpOptions->padfDstNoDataImag = (double*) CPLMalloc(channels() * sizeof(double));
for (unsigned int i = 0; i < channels(); ++i) {
int hasNodataVal;
double nodataVal = poSrcDS->GetRasterBand(i+1)->GetNoDataValue(&hasNodataVal);
psWarpOptions->padfSrcNoDataReal[i] = nodataVal;
psWarpOptions->padfSrcNoDataImag[i] = 0;
if (to.hasNodataValue())
nodataVal = poDstDS->GetRasterBand(i+1)->GetNoDataValue(&hasNodataVal);
psWarpOptions->padfDstNoDataReal[i] = nodataVal;
psWarpOptions->padfDstNoDataImag[i] = 0;
}
}
psWarpOptions->papszWarpOptions = nullptr;
if (channels() == 1)
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "UNIFIED_SRC_NODATA", "YES" ); // otherwise no data points will be replaced by the nearest valid neighbor
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "SAMPLE_STEPS", "32"); // increase precision to some 2^x number
// workaround for a bug in GDAL, which sets the scale factor for an identity warp to 0.5 instead of 1 and then uses anti-aliasing
if (from.geotransSRS.IsSame(&to.geotransSRS) &&
from.geotrans.XToX == to.geotrans.XToX &&
from.geotrans.YToY == to.geotrans.YToY &&
from.geotrans.XToY == to.geotrans.XToY &&
from.geotrans.YToX == to.geotrans.YToX)
{
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "XSCALE", "1");
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "YSCALE", "1");
}
// Establish reprojection transformer.
char** papszTO = nullptr;
papszTO = CSLSetNameValue(papszTO, "STRIP_VERT_CS", "YES");
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer2(poSrcDS, poDstDS, papszTO);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// Initialize and execute the warp operation.
GDALWarpOperation oOperation;
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0, 0, sz.width, sz.height);
GDALClose(poDstDS);
GDALClose(poSrcDS);
// workaround for a bug in GDAL, which uses the nearest valid neighbor, even if the nearest neighbor is invalid. This happens only if not using nearest neighbor interpolation and UNIFIED_SRC_NODATA is set to NO, i. e. for multi-channel images
if (channels() > 1 && method != InterpMethod::nearest && from.hasNodataValue()) {
// make mask from source nodata values
std::vector<Interval> nodataValIntervals;
std::vector<double> nodataVals;
for (unsigned int i = 0; i < channels(); ++i) {
double d = from.getNodataValue(i);
nodataVals.push_back(d);
nodataValIntervals.push_back(Interval::closed(d, d));
}
Image mask = createMultiChannelMaskFromRange(nodataValIntervals);
// warp source nodata mask with nearest neighbor interpolation
Image warpedMask{sz, mask.type()};
GDALDataset* poSrcMaskDS = const_cast<GDALDataset*>(mask.asGDALDataset());
GDALDataset* poDstMaskDS = warpedMask.asGDALDataset();
from.addTo(poSrcMaskDS);
to.addTo(poDstMaskDS);
psWarpOptions->eResampleAlg = GRA_NearestNeighbour;
psWarpOptions->hSrcDS = poSrcMaskDS;
psWarpOptions->hDstDS = poDstMaskDS;
for (unsigned int i = 0; i < channels(); ++i) { psWarpOptions->padfSrcNoDataReal[i] = std::numeric_limits<double>::quiet_NaN();;
psWarpOptions->padfDstNoDataReal[i] = std::numeric_limits<double>::quiet_NaN();;
}
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0, 0, sz.width, sz.height);
GDALClose(poDstMaskDS);
GDALClose(poSrcMaskDS);
// set dst nodata values
if (to.hasNodataValue())
for (unsigned int i = 0; i < channels(); ++i)
nodataVals.at(i) = to.getNodataValue(i);
warped.set(nodataVals, warpedMask);
}
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
GDALDestroyWarpOptions(psWarpOptions);
return warped;
}
std::vector<Image> ConstImage::split(std::vector<unsigned int> chans) const {
std::vector<cv::Mat> cvImages;
if (chans.empty()) {
// simply split all channels
cvImages = std::vector<cv::Mat>(channels());
cv::split(img, cvImages);
}
else {
// from channel - to channel relation
std::vector<int> fromTo;
for (unsigned int i = 0; i < chans.size(); ++i) {
unsigned int ch = chans[i];
if (ch >= channels())
IF_THROW_EXCEPTION(image_type_error("You requested to split an image and to use channel "
+ std::to_string(ch) + ", but the image has got only "
+ std::to_string(channels()))) << errinfo_image_type(type());
fromTo.push_back(ch);
fromTo.push_back(i);
}
// prepare single channel output images
cvImages = std::vector<cv::Mat>(chans.size());
for (cv::Mat& m : cvImages)
m.create(size(), toCVType(basetype()));
// split
std::vector<cv::Mat> src{cvMat()};
cv::mixChannels(src, cvImages, fromTo.data(), chans.size());
}
// convert to Image type
std::vector<Image> channelImages(cvImages.size());
for (unsigned int i = 0; i < channelImages.size(); ++i)
channelImages[i].img = std::move(cvImages[i]);
return channelImages;
}
Image merge(std::vector<Image> const& images) {
return merge_impl(images);
}
Image merge(std::vector<ConstImage> const& images) {
return merge_impl(images);
}
namespace {template<typename OP>
Image SimpleBroadcastOperation(ConstImage const& A, Image B, OP const& op) {
cv::Mat a = A.cvMat();
if (a.channels() == 1 && B.channels() > 1) {
std::vector<cv::Mat> aaa(B.channels(), a);
cv::merge(aaa, a); // note: broadcasting could be done without temporary memory allocation for better performance
}
else if (B.channels() == 1 && a.channels() > 1) {
std::vector<cv::Mat> bbb(a.channels(), B.cvMat());
cv::merge(bbb, B.cvMat());
}
op(a, B.cvMat(), B.cvMat());
return B;
}
} /* anonymous namespace */
Image ConstImage::absdiff(Image&& B) const& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::absdiff(mat1, mat2, res);});
}
Image ConstImage::absdiff(ConstImage const& B) const& {
if (B.channels() > channels()) // for performance only in case of broadcasting
return absdiff(B.clone()); // B has multiple channels and can be used to save the result
return B.absdiff(clone());
}
Image Image::absdiff(ConstImage const& B)&& {
return B.absdiff(std::move(*this));
}
Image ConstImage::abs() const& {
if (isUnsignedType(type()))
return clone();
if (channels() <= 4) {
Image ret;
ret.img = cv::abs(img);
return ret;
}
Image ret = clone();
auto op = [] (auto& v, int, int, int) {
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(std::abs(v));
};
CallBaseTypeFunctorRestrictBaseTypesTo<Type::int8, Type::int16, Type::int32, Type::float32, Type::float64>::run(InplacePointOperationFunctor{ret, {}, op}, ret.basetype());
return ret;
}
Image Image::abs()&& {
if (isUnsignedType(type()))
return std::move(*this);
if (channels() <= 4) {
img = cv::abs(img);
}
else {
auto op = [] (auto& v, int, int, int) {
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(std::abs(v));
};
CallBaseTypeFunctorRestrictBaseTypesTo<Type::int8, Type::int16, Type::int32, Type::float32, Type::float64>::run(InplacePointOperationFunctor{*this, {}, op}, basetype()); }
return std::move(*this);
}
namespace {
void checkMask(ConstImage const& mask, unsigned int imgChans) {
if (mask.channels() > 1 && mask.channels() != imgChans)
IF_THROW_EXCEPTION(image_type_error("Mask has a wrong number of channels. It should be single-channel"
+ (imgChans > 1 ? " or multi-channel with " + std::to_string(imgChans) + " channels" : "")
+ ". However, it has " + std::to_string(mask.channels()) + " channels."))
<< errinfo_image_type(mask.type());
if (mask.basetype() != Type::uint8)
IF_THROW_EXCEPTION(image_type_error("Mask has an invalid type. It should have basetype Type::uint8, "
"but has " + to_string(mask.basetype()) + "."))
<< errinfo_image_type(mask.basetype());
}
template<typename CvOp, typename StdOp>
Image minimum_maximum(Image A, ConstImage const& B, ConstImage const& mask, CvOp const& cvminmax, StdOp const& stdminmax) {
if (B.empty())
return A;
if (A.empty())
return B.clone();
// no mask --> use OpenCV with broadcasting
if (mask.empty())
return SimpleBroadcastOperation(B, std::move(A), cvminmax);
// broadcasting minimum or maximum with mask
if (A.channels() > 1 && B.channels() == 1) {
auto op = [&B, &stdminmax](auto& v, int const& x, int const& y, int const& /*c*/){
using type = std::remove_reference_t<decltype(v)>;
v = stdminmax(v, B.at<type>(x, y, 0));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
return A;
}
if (A.channels() == 1 && B.channels() > 1) {
std::vector<cv::Mat> aaa(B.channels(), A.cvMat());
cv::merge(aaa, A.cvMat());
}
auto op = [&B, &stdminmax](auto& v, int const& x, int const& y, int const& c){
using type = std::remove_reference_t<decltype(v)>;
v = stdminmax(v, B.at<type>(x, y, c));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
return A;
}
template<typename CvOp, typename CustomOp>
Image minimum_maximum_pixel(std::vector<double> const& pix, Image A, ConstImage const& mask, CvOp const& cvop, CustomOp const& customop) {
using std::to_string;
checkMask(mask, A.channels());
if (pix.size() != 1 && pix.size() != A.channels())
IF_THROW_EXCEPTION(invalid_argument_error("The number of pixel values must be one or match number of channels of the image. "
"Here the number of values is " + to_string(pix.size()) + ", but the image has " + to_string(A.channels()) + "."));
// no mask --> use OpenCV
if (mask.empty()) {
cvop(A.cvMat(), pix, A.cvMat());
return A;
}
// fallback code
if (pix.size() < A.channels()) {
double p = pix.front();
auto op = [p, &customop] (auto& v, int /*x*/, int /*y*/, int /*c*/) {
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(customop(static_cast<double>(v), p));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
}
else {
auto op = [&pix, &customop](auto& v, int, int, int const& c){
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(customop(static_cast<double>(v), pix.at(c)));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
}
return A;
}
template<typename CvOp, typename CustomOp>
Image add_subtract_pixel(std::vector<double> const& pix, Image A, ConstImage const& mask, CvOp const& cvop, CustomOp const& customop) {
if (A.empty())
return A;
using std::to_string;
checkMask(mask, A.channels());
if (pix.size() != 1 && pix.size() != A.channels())
IF_THROW_EXCEPTION(invalid_argument_error("The number of pixel values must be one or match number of channels of the image. "
"Here the number of values is " + to_string(pix.size()) + ", but the image has " + to_string(A.channels()) + "."));
// try OpenCV operation first, for performance reasons, but it does not accept all types (e. g. it throws on uint8x2)
if (mask.empty() || mask.channels() == 1) {
try {
cvop(A.cvMat(), pix, A.cvMat(), mask.cvMat());
return A;
} catch (cv::Exception&) {
/* empty, fall through */
}
}
// fallback code
if (pix.size() < A.channels()) {
double p = pix.front();
auto op = [p, &customop] (auto& v, int /*x*/, int /*y*/, int /*c*/) {
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(customop(v, p));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
}
else {
auto op = [&pix, &customop](auto& v, int, int, int const& c){
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(customop(v, pix.at(c)));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{A, mask, op}, A.type());
}
return A;
}
} /* anonymous namespace */
Image ConstImage::minimum(Image&& B, ConstImage const& mask) const& {
if (mask.empty())
return std::move(B).minimum(*this, mask);
return clone().minimum(B, mask); // order must be consistent when using masks
}
Image ConstImage::minimum(ConstImage const& B, ConstImage const& mask) const& {
if (B.channels() < channels() || !mask.empty()) // for performance only in case of broadcasting and non-empty mask
return clone().minimum(B, mask); // A has multiple channels and can be used to save the result
return B.clone().minimum(*this, mask); // in case of a non-empty mask it must be A.minimum(B, mask) instead of B.minimum(A, mask)
}
Image Image::minimum(ConstImage const& B, ConstImage const& mask)&& {
return minimum_maximum(std::move(*this), B, mask,
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::min(mat1, mat2, res);},
[] (auto a, auto b) { return std::min(a, b); });
}
Image ConstImage::minimum(std::vector<double> const& pix, ConstImage const& mask) const& {
return clone().minimum(pix, mask); // call rvalue method
}
Image Image::minimum(std::vector<double> const& pix, ConstImage const& mask)&& {
return minimum_maximum_pixel(pix, std::move(*this), mask,
[] (cv::Mat const& mat, std::vector<double> const& pix, cv::Mat& res) {
cv::min(mat, pix, res);
},
[] (auto a, auto b) { return std::min(a, b); });
}
Image ConstImage::maximum(Image&& B, ConstImage const& mask) const& {
if (mask.empty())
return std::move(B).maximum(*this, mask);
return clone().maximum(B, mask); // order must be consistent when using masks
}
Image ConstImage::maximum(ConstImage const& B, ConstImage const& mask) const& {
if (B.channels() < channels() || !mask.empty()) // for performance only in case of broadcasting and non-empty mask
return clone().maximum(B, mask); // A has multiple channels and can be used to save the result
return B.clone().maximum(*this, mask); // in case of a non-empty mask it must be A.maximum(B, mask) instead of B.maximum(A, mask)
}
Image Image::maximum(ConstImage const& B, ConstImage const& mask)&& {
return minimum_maximum(std::move(*this), B, mask,
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::max(mat1, mat2, res);},
[] (auto a, auto b) { return std::max(a, b); });
}
Image ConstImage::maximum(std::vector<double> const& pix, ConstImage const& mask) const& {
return clone().maximum(pix, mask); // call rvalue method
}
Image Image::maximum(std::vector<double> const& pix, ConstImage const& mask)&& {
return minimum_maximum_pixel(pix, std::move(*this), mask,
[] (cv::Mat const& mat, std::vector<double> const& pix, cv::Mat& res) {
cv::max(mat, pix, res);
},
[] (auto a, auto b) { return std::max(a, b); });
}
Image ConstImage::add(ConstImage const& B, Type resultType) const {
Image ret;
cv::add(img, B.img, ret.img, cv::noArray(), toCVType(resultType));
return ret;
}
Image ConstImage::add(Image&& B) const& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::add(mat1, mat2, res);});
}
Image ConstImage::add(ConstImage const& B) const& {
if (B.channels() > channels()) // for performance only in case of broadcasting
return add(B.clone()); // B has multiple channels and can be used to save the result
return B.add(clone());
}
Image Image::add(ConstImage const& B)&& {
return B.add(std::move(*this));
}
Image ConstImage::add(std::vector<double> const& pix, ConstImage const& mask, Type resultType) const& {
if (resultType == Type::invalid || getBaseType(resultType) == basetype())
return clone().add(pix, mask); // call rvalue method
else // different type
return convertTo(resultType).add(pix, mask); // call rvalue method
}
Image Image::add(std::vector<double> const& pix, ConstImage const& mask)&& {
int sameResultType = toCVType(basetype());
return add_subtract_pixel(pix, std::move(*this), mask,
[sameResultType] (cv::Mat const& mat, std::vector<double> const& pix, cv::Mat& res, cv::Mat const& mask) {
cv::add(mat, pix, res, mask, sameResultType);
},
[] (auto a, auto b) { return a + b; });
}
Image ConstImage::subtract(ConstImage const& B, Type resultType) const {
Image ret;
cv::subtract(img, B.img, ret.img, cv::noArray(), toCVType(resultType));
return ret;
}
Image ConstImage::subtract(Image&& B) const& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(*this, std::move(B),
[](cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res){cv::subtract(mat1, mat2, res);});
}
Image ConstImage::subtract(ConstImage const& B) const& {
return subtract(B.clone());
}
Image Image::subtract(ConstImage const& B)&& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(/*mat1*/ B, /*mat2*/ std::move(*this),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::subtract(mat2, mat1, res);});
}
Image ConstImage::subtract(std::vector<double> const& pix, ConstImage const& mask, Type resultType) const& {
if (resultType == Type::invalid || getBaseType(resultType) == basetype())
return clone().subtract(pix, mask); // call rvalue method
else // different type
return convertTo(resultType).subtract(pix, mask); // call rvalue method
}
Image Image::subtract(std::vector<double> const& pix, ConstImage const& mask)&& {
int sameResultType = toCVType(basetype());
return add_subtract_pixel(pix, std::move(*this), mask,
[sameResultType] (cv::Mat const& mat, std::vector<double> const& pix, cv::Mat& res, cv::Mat const& mask) {
cv::subtract(mat, pix, res, mask, sameResultType);
},
[] (auto a, auto b) { return a - b; });
}
Image ConstImage::multiply(ConstImage const& B, Type resultType) const {
Image ret;
cv::multiply(img, B.img, ret.img, 1, toCVType(resultType));
return ret;
}
Image ConstImage::multiply(Image&& B) const& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::multiply(mat1, mat2, res);});
}
Image ConstImage::multiply(ConstImage const& B) const& {
if (B.channels() > channels()) // for performance only in case of broadcasting
return multiply(B.clone()); // B has multiple channels and can be used to save the result
return B.multiply(clone());
}
Image Image::multiply(ConstImage const& B)&& {
return B.multiply(std::move(*this));
}
Image ConstImage::multiply(std::vector<double> const& pix, ConstImage const& mask, Type resultType) const& {
if (resultType == Type::invalid || getBaseType(resultType) == basetype())
return clone().multiply(pix, mask); // call rvalue method
else // different type
return convertTo(resultType).multiply(pix, mask); // call rvalue method
}
Image Image::multiply(std::vector<double> const& pix, ConstImage const& mask)&& {
if (empty())
return std::move(*this);
using std::to_string;
checkMask(mask, channels());
if (pix.size() != 1 && pix.size() != channels())
IF_THROW_EXCEPTION(invalid_argument_error("The number of pixel values must be one or match number of channels of the image. "
"Here the number of values is " + to_string(pix.size()) + ", but the image has " + to_string(channels()) + "."));
// try OpenCV operation first, for performance reasons, but it does not accept all types (e. g. it throws on uint8x2)
if (mask.empty()) {
try {
cv::multiply(img, pix, img, /* scale */ 1, toCVType(basetype()));
return std::move(*this);
} catch (cv::Exception&) {
/* empty, fall through */
}
}
std::vector<double> const& pix_n = pix.size() < channels() ? std::vector<double>(channels(), pix.front()) : pix;
auto op = [&pix_n](auto& v, int, int, int const& c){
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(v * pix_n.at(c));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{*this, mask, op}, basetype());
return std::move(*this);
}
Image ConstImage::divide(ConstImage const& B, Type resultType) const {
Image ret;
cv::divide(img, B.img, ret.img, 1, toCVType(resultType));
return ret;
}
Image ConstImage::divide(Image&& B) const& { if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::divide(mat1, mat2, res);});
}
Image ConstImage::divide(ConstImage const& B) const& {
return divide(B.clone());
}
Image Image::divide(ConstImage const& B)&& {
if (empty() || B.empty())
IF_THROW_EXCEPTION(invalid_argument_error("Empty images are not allowed for arithmetic operations (only for bitwise logical operations)."));
return SimpleBroadcastOperation(/*mat1*/ B, /*mat2*/ std::move(*this),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::divide(mat2, mat1, res);});
}
Image ConstImage::divide(std::vector<double> const& pix, ConstImage const& mask, Type resultType) const& {
if (resultType == Type::invalid || getBaseType(resultType) == basetype())
return clone().divide(pix, mask); // call rvalue method
else // different type
return convertTo(resultType).divide(pix, mask); // call rvalue method
}
Image Image::divide(std::vector<double> const& pix, ConstImage const& mask)&& {
if (empty())
return std::move(*this);
using std::to_string;
checkMask(mask, channels());
if (pix.size() != 1 && pix.size() != channels())
IF_THROW_EXCEPTION(invalid_argument_error("The number of pixel values must be one or match number of channels of the image. "
"Here the number of values is " + to_string(pix.size()) + ", but the image has " + to_string(channels()) + "."));
for (double d : pix) {
if (d == 0) {
std::string vals = "[";
for (unsigned int i = 0; i < pix.size(); ++i) {
vals += std::to_string(pix.at(i));
if (i < pix.size() - 1)
vals += ", ";
else
vals += "]";
}
IF_THROW_EXCEPTION(invalid_argument_error("You try to divide the image by zero. Your divisor is: " + vals));
}
}
// try OpenCV operation first, for performance reasons, but it does not accept all types (e. g. it throws on uint8x2)
if (mask.empty()) {
try {
cv::divide(img, pix, img, /* scale */ 1, toCVType(basetype()));
return std::move(*this);
} catch (cv::Exception&) {
/* empty, fall through */
}
}
std::vector<double> const& pix_n = pix.size() < channels() ? std::vector<double>(channels(), pix.front()) : pix;
auto op = [&pix_n](auto& v, int, int, int const& c){
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(v / pix_n.at(c));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{*this, mask, op}, type());
return std::move(*this);}
Image ConstImage::bitwise_and(Image&& B) const& {
if (B.empty())
return clone();
if (empty())
return std::move(B);
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::bitwise_and(mat1, mat2, res);});
}
Image ConstImage::bitwise_and(ConstImage const& B) const& {
return bitwise_and(B.clone());
}
Image Image::bitwise_and(ConstImage const& B)&& {
return B.bitwise_and(std::move(*this));
}
Image ConstImage::bitwise_or(Image&& B) const& {
if (B.empty())
return clone();
if (empty())
return std::move(B);
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::bitwise_or(mat1, mat2, res);});
}
Image ConstImage::bitwise_or(ConstImage const& B) const& {
return bitwise_or(B.clone());
}
Image Image::bitwise_or(ConstImage const& B)&& {
return B.bitwise_or(std::move(*this));
}
Image ConstImage::bitwise_xor(Image&& B) const& {
if (B.empty())
return clone();
if (empty())
return std::move(B);
return SimpleBroadcastOperation(*this, std::move(B),
[] (cv::Mat const& mat1, cv::Mat const& mat2, cv::Mat& res) {cv::bitwise_xor(mat1, mat2, res);});
}
Image ConstImage::bitwise_xor(ConstImage const& B) const& {
return bitwise_xor(B.clone());
}
Image Image::bitwise_xor(ConstImage const& B)&& {
return B.bitwise_xor(std::move(*this));
}
Image ConstImage::bitwise_not() const& {
Image ret;
if (!empty())
cv::bitwise_not(img, ret.img);
return ret;
}
Image Image::bitwise_not()&& {
if (!empty())
cv::bitwise_not(img, img); return std::move(*this);
}
std::vector<std::pair<ValueWithLocation, ValueWithLocation>> ConstImage::minMaxLocations(ConstImage const& mask) const {
std::vector<std::pair<ValueWithLocation, ValueWithLocation>> minmax(channels());
if(mask.channels() != 1 && mask.channels() != channels())
IF_THROW_EXCEPTION(image_type_error("The mask must have one channel or the same number of channels as the image. The image has "
+ std::to_string(channels()) + " and the mask " + std::to_string(mask.channels()) + " channels."))
<< errinfo_image_type(mask.type());
auto img_layers = split();
auto mask_layers = mask.split();
for (unsigned int c = 0; c < channels(); ++c) {
cv::Mat const& img = img_layers.at(c).img;
double min, max;
cv::Point pMin, pMax;
if (mask.empty())
cv::minMaxLoc(img, &min, &max, &pMin, &pMax);
else if (mask.channels() == 1) {
cv::minMaxLoc(img, &min, &max, &pMin, &pMax, mask.img);
}
else { // (mask.channels() == channels)
cv::Mat const& m = mask_layers.at(c).img;
cv::minMaxLoc(img, &min, &max, &pMin, &pMax, m);
}
minmax.at(c) = std::make_pair(ValueWithLocation{min, pMin}, ValueWithLocation{max, pMax});
}
return minmax;
}
std::vector<double> ConstImage::mean(ConstImage const& mask) const {
std::vector<double> mean(channels());
if (mask.channels() == 1 && channels() <= 4) {
cv::Scalar mean_cv = cv::mean(img, mask.img);
for (unsigned int i = 0; i < channels(); ++i)
mean.at(i) = mean_cv[i];
}
else {
if(mask.channels() != 1 && mask.channels() != channels())
IF_THROW_EXCEPTION(image_type_error("The mask must have one channel or the same number of channels as the image"))
<< errinfo_image_type(mask.type());
auto img_layers = split();
auto mask_layers = mask.split();
cv::Scalar temp_mean;
for (unsigned int c = 0; c < channels(); ++c) {
cv::Mat const& img = img_layers.at(c).img;
if (mask.empty())
temp_mean = cv::mean(img);
else if (mask.channels() == 1)
temp_mean = cv::mean(img, mask.img);
else { // (mask.channels() == channels)
cv::Mat const& m = mask_layers.at(c).img;
temp_mean = cv::mean(img, m);
}
mean.at(c) = temp_mean[0];
}
}
return mean;
}
std::pair<std::vector<double>, std::vector<double>> ConstImage::meanStdDev(ConstImage const& mask, bool sampleCorrection) const {
std::vector<double> mean(channels());
std::vector<double> stdDev(channels());
if (mask.channels() == 1 && channels() <= 4) {
double cor = 1; if (sampleCorrection) {
double n = mask.empty() ? size().area() : cv::countNonZero(mask.img);
if (n > 1)
cor = std::sqrt(n / (n - 1));
}
cv::Scalar mean_cv, stddev_cv;
cv::meanStdDev(img, mean_cv, stddev_cv, mask.img);
for (unsigned int i = 0; i < channels(); ++i) {
mean.at(i) = mean_cv[i];
stdDev.at(i) = stddev_cv[i] * cor;
}
}
else {
if(mask.channels() != 1 && mask.channels() != channels())
IF_THROW_EXCEPTION(image_type_error("The mask must have one channel or the same number of channels as the image"))
<< errinfo_image_type(mask.type());
auto img_layers = split();
auto mask_layers = mask.split();
cv::Scalar temp_mean, temp_stddev;
for (unsigned int c = 0; c < channels(); ++c) {
cv::Mat const& img = img_layers.at(c).img;
double n = size().area();
if (mask.empty())
cv::meanStdDev(img, temp_mean, temp_stddev);
else if (mask.channels() == 1) {
cv::meanStdDev(img, temp_mean, temp_stddev, mask.img);
if (sampleCorrection)
n = cv::countNonZero(mask.img);
}
else { // (mask.channels() == channels)
cv::Mat const& m = mask_layers.at(c).img;
cv::meanStdDev(img, temp_mean, temp_stddev, m);
if (sampleCorrection)
n = cv::countNonZero(m);
}
double cor = 1;
if (sampleCorrection && n > 1)
cor = std::sqrt(n / (n - 1));
mean.at(c) = temp_mean[0];
stdDev.at(c) = temp_stddev[0] * cor;
}
}
return {mean, stdDev};
}
namespace {
template<class ForwardImgIt>
ForwardImgIt unique_with_mask(ForwardImgIt img_first, ForwardImgIt img_last, cv::MatConstIterator_<uint8_t> mask_first, cv::MatConstIterator_<uint8_t> mask_last)
{
if (img_first == img_last || mask_first == mask_last)
return img_last;
ForwardImgIt img_result = img_first;
while (!*mask_first && mask_first != mask_last) {
++mask_first;
++img_first;
}
// check if no location was valid
if (mask_first == mask_last)
return img_result;
// check if first location was invalid
if (img_result != img_first)
// make sure the first position is valid
*img_result = std::move(*img_first);
// iterate
while (++img_first != img_last && ++mask_first != mask_last)
if (*mask_first && !(*img_result == *img_first) && ++img_result != img_first)
*img_result = std::move(*img_first);
return ++img_result;
}
struct UniqueFunctor {
Image& img;
ConstImage const& mask;
template<Type basetype>
std::vector<double> operator()() const {
static_assert(getChannels(basetype) == 1, "This functor only accepts base type to reduce code size.");
if (img.channels() != 1)
IF_THROW_EXCEPTION(invalid_argument_error("Only single-channel images supported for unique. The image has " + std::to_string(img.channels()) + " channels.")) << errinfo_image_type(img.type());
if (mask.type() != Type::uint8x1)
IF_THROW_EXCEPTION(invalid_argument_error("Only single-channel masks (with Type uint8) supported for unique. The mask has type " + to_string(mask.type()) + ".")) << errinfo_image_type(mask.type());
if (!mask.empty() && mask.size() != img.size())
IF_THROW_EXCEPTION(invalid_argument_error("Mask and image have different sizes. image: " + to_string(img.size()) + ", mask: " + to_string(mask.size())));
using type = typename DataType<basetype>::base_type;
auto begin = img.begin<type>(), end = img.end<type>();
// remove adjacent duplicates to reduce size
auto last = mask.empty() ? std::unique(begin, end) : unique_with_mask(begin, end, mask.begin<uint8_t>(), mask.end<uint8_t>());
std::sort(begin, last); // sort remaining elements
last = std::unique(begin, last); // remove duplicates
std::vector<double> vals;
while (begin != last) {
vals.push_back(static_cast<double>(*begin));
++begin;
}
return vals;
}
};
struct CountUniqueFunctor {
ConstImage const& img;
ConstImage const& mask;
template<Type basetype>
std::map<double, unsigned int> operator()() const {
static_assert(getChannels(basetype) == 1, "This functor only accepts base type to reduce code size.");
if (img.channels() != 1)
IF_THROW_EXCEPTION(invalid_argument_error("Only single-channel images supported for unique. The image has " + std::to_string(img.channels()) + " channels.")) << errinfo_image_type(img.type());
if (mask.type() != Type::uint8x1)
IF_THROW_EXCEPTION(invalid_argument_error("Only single-channel masks (with Type uint8) supported for unique. The mask has type " + to_string(mask.type()) + ".")) << errinfo_image_type(mask.type());
if (!mask.empty() && mask.size() != img.size())
IF_THROW_EXCEPTION(invalid_argument_error("Mask and image have different sizes. image: " + to_string(img.size()) + ", mask: " + to_string(mask.size())));
using type = typename DataType<basetype>::base_type;
std::map<double, unsigned int> unique;
auto img_it = img.begin<type>(), img_end = img.end<type>();
if (mask.empty()) {
for (; img_it != img_end; ++img_it)
++unique[*img_it];
}
else {
auto mask_it = mask.begin<uint8_t>(), mask_end = mask.end<uint8_t>();
for (; img_it != img_end && mask_it != mask_end; ++img_it, ++mask_it)
if (*mask_it)
++unique[*img_it];
} return unique;
}
};
} /* anonymous namespace */
std::vector<double> ConstImage::unique(ConstImage const& mask) const& {
return clone().unique(mask);
}
std::vector<double> Image::unique(ConstImage const& mask)&& {
return CallBaseTypeFunctor::run(UniqueFunctor{*this, mask}, basetype());
}
std::map<double, unsigned int> ConstImage::uniqueWithCount(ConstImage const& mask) const {
return CallBaseTypeFunctor::run(CountUniqueFunctor{*this, mask}, basetype());
}
Image& Image::set(double val, ConstImage const& mask) {
return set(std::vector<double>(channels(), val), mask);
}
Image& Image::set(std::vector<double> const& vals, ConstImage const& mask) {
using std::to_string;
checkMask(mask, channels());
if (vals.size() != channels())
IF_THROW_EXCEPTION(invalid_argument_error("The number of values must match number of channels as the image. "
"Here the number of values is " + to_string(vals.size()) + ", but the image has " + to_string(channels()) + "."));
if (mask.type() == Type::uint8x1 && channels() <= 4) { // OpenCV supports multi-channel masks since version 3.3.1. So in future, we can simplify this.
cv::Scalar cvVals;
for (unsigned int i = 0; i < vals.size() && i < 4; ++i)
cvVals[i] = vals.at(i);
img.setTo(cvVals, mask.cvMat());
return *this;
}
auto op = [&vals] (auto& v, int, int, int const& c) {
using type = std::remove_reference_t<decltype(v)>;
v = cv::saturate_cast<type>(vals.at(c));
};
CallBaseTypeFunctor::run(InplacePointOperationFunctor{*this, mask, op}, basetype());
// Alternative implementation based on OpenCV functions. Performance results for an older version:
// Requires more memory (ca. x2), could be faster for some situations and slower for other (33% faster or slower).
// unsigned int chans = channels();
// std::vector<cv::Mat> img_chans(chans);
// cv::split(img, img_chans);
// std::vector<cv::Mat> mask_chans(chans);
// cv::split(mask.img, mask_chans);
// for (unsigned int c = 0; c < chans; ++c)
// img_chans[c].setTo(val[c], mask_chans[c]);
// cv::merge(img_chans, img);
return *this;
}
namespace {
double tolerance(double v, Type img_type) {
if (getBaseType(img_type) == Type::float32)
return 2 * std::numeric_limits<float>::epsilon() * std::abs(v);
else // (getBaseType(img_type) == Type::float64)
return 2 * std::numeric_limits<double>::epsilon() * std::abs(v);
}
}
static std::array<std::vector<double>,2> fixBounds(std::vector<Interval> const& channelRanges, Type t) {
using std::to_string;
int size = channelRanges.size();
int chans = getChannels(t);
if (size != 1 && size != chans)
IF_THROW_EXCEPTION(image_type_error("The number of ranges is not valid. Either give one range that is used for all channels or give one range per channel. "
"Here the image has " + to_string(chans) + " channels and the number of ranges is " + to_string(size) + "."))
<< errinfo_image_type(t);
// for integer images, values must be in the int32 range
std::vector<double> lowBounds;
std::vector<double> highBounds;
if (isIntegerType(t)) {
for (Interval const& i : channelRanges) {
// limit on int32 range, but keep nan
bool leftClosed = i.bounds().left() != boost::icl::interval_bounds::open();
bool rightClosed = i.bounds().right() != boost::icl::interval_bounds::open();
double l = i.lower();
double u = i.upper();
if (l < std::numeric_limits<int32_t>::min()) {
l = std::numeric_limits<int32_t>::min();
leftClosed = true;
}
if (u > std::numeric_limits<int32_t>::max()) {
u = std::numeric_limits<int32_t>::max();
rightClosed = true;
}
// handle open intervals
if (leftClosed)
l = std::ceil(l);
else
l = std::floor(l);
if (rightClosed)
u = std::floor(u);
else
u = std::ceil(u);
lowBounds.push_back( std::isnan(l) || leftClosed ? l : l + 1);
highBounds.push_back(std::isnan(u) || rightClosed ? u : u - 1);
}
}
else { // floating point image. Handle open intervals by using tolerances with 2 ULP (units in the last place)
for (Interval const& i : channelRanges) {
double l = i.lower();
double u = i.upper();
if (i.bounds().left() == boost::icl::interval_bounds::open())
l += tolerance(l, t);
if (i.bounds().right() == boost::icl::interval_bounds::open())
u -= tolerance(u, t);
lowBounds.push_back(l);
highBounds.push_back(u);
}
}
// fill up bounds
if (size < chans) {
lowBounds.resize(chans, lowBounds.front());
highBounds.resize(chans, highBounds.front());
}
return {std::move(lowBounds), std::move(highBounds)};
}
static std::vector<IntervalSet> fixBounds(std::vector<IntervalSet> channelSets, Type t) {
using std::to_string;
int size = channelSets.size();
int chans = getChannels(t);
if (size != 1 && size != chans) IF_THROW_EXCEPTION(image_type_error("The number of sets is not valid. Either give one set that is used for all channels or give one set per channel. "
"Here the image has " + to_string(chans) + " channels and the number of sets is " + to_string(size) + "."))
<< errinfo_image_type(t);
// for integer images, values must be in the int32 range and we like to convert intervals to closed intervals
if (isIntegerType(t))
for (IntervalSet& is : channelSets)
is = discretizeBounds(is);
else { // floating point image. Handle open intervals by using tolerances with 2 ULP (units in the last place)
for (IntervalSet& is : channelSets) {
IntervalSet interval_set_with_tolerances;
for (Interval const& i : is) {
double l = i.lower();
double u = i.upper();
if (i.bounds().left() == boost::icl::interval_bounds::open())
l += tolerance(l, t);
if (i.bounds().right() == boost::icl::interval_bounds::open())
u -= tolerance(u, t);
interval_set_with_tolerances += Interval(l, u, i.bounds());
}
is = interval_set_with_tolerances;
}
}
// fill up sets
if (size < chans)
channelSets.resize(chans, channelSets.front());
return channelSets;
}
Image ConstImage::createSingleChannelMaskFromRange(std::vector<Interval> const& channelRanges, bool useAnd) const {
auto l_h = fixBounds(channelRanges, type());
std::vector<double> const& low = l_h[0];
std::vector<double> const& high = l_h[1];
Image ret;
if (useAnd)
cv::inRange(img, low, high, ret.img);
else {
// apply inRange channel-wise
unsigned int chans = channels();
std::vector<cv::Mat> img_chans(chans);
cv::split(img, img_chans);
for (unsigned int c = 0; c < chans; ++c) {
if (c == 0)
// first iteration
cv::inRange(img_chans[c], low[c], high[c], ret.cvMat());
else {
cv::Mat temp_mask;
cv::inRange(img_chans[c], low[c], high[c], temp_mask);
ret.img |= temp_mask;
}
}
}
return ret;
}
Image ConstImage::createSingleChannelMaskFromRange(Interval const& channelRange, bool useAnd) const {
return createSingleChannelMaskFromRange(std::vector<Interval>{channelRange}, useAnd);
}
Image ConstImage::createMultiChannelMaskFromRange(std::vector<Interval> const& channelRanges) const {
unsigned int chans = channels();
if (chans == 1)
return createSingleChannelMaskFromRange(channelRanges);
auto l_h = fixBounds(channelRanges, type());
std::vector<double> const& low = l_h[0]; std::vector<double> const& high = l_h[1];
// apply inRange channel-wise
std::vector<cv::Mat> img_chans(chans);
cv::split(img, img_chans);
std::vector<cv::Mat> mask_chans(chans);
for (unsigned int c = 0; c < chans; ++c)
cv::inRange(img_chans[c], low[c], high[c], mask_chans[c]);
Image ret;
cv::merge(mask_chans, ret.img);
return ret;
}
Image ConstImage::createMultiChannelMaskFromRange(Interval const& channelRange) const {
return createMultiChannelMaskFromRange(std::vector<Interval>{channelRange});
}
static std::vector<cv::Mat> getMaskVec(ConstImage const& src, std::vector<IntervalSet> const& channelSets) {
assert(static_cast<unsigned int>(src.channels()) == channelSets.size());
// apply inRange channel- and interval-wise, use OR to relate the intervals
unsigned int chans = src.channels();
std::vector<cv::Mat> img_chans(chans);
cv::split(src.cvMat(), img_chans);
std::vector<cv::Mat> mask_chans(chans);
cv::Mat temp;
for (unsigned int c = 0; c < chans; ++c) {
mask_chans[c] = cv::Mat(src.height(), src.width(), CV_8UC1);
mask_chans[c] = cv::Scalar::all(0);
IntervalSet const& is = channelSets.at(c);
for (Interval const& i : is) {
cv::inRange(img_chans[c], i.lower(), i.upper(), temp);
mask_chans[c] |= temp;
}
}
return mask_chans;
}
Image ConstImage::createSingleChannelMaskFromSet(std::vector<IntervalSet> const& channelSets, bool useAnd) const {
auto sets = fixBounds(channelSets, type());
std::vector<cv::Mat> mask_chans = getMaskVec(*this, sets);
// combine channel masks with AND
Image ret(size(), Type::uint8x1);
if (useAnd) {
ret.set(255);
for (cv::Mat const& m : mask_chans)
ret.img &= m;
}
else {
ret.set(0);
for (cv::Mat const& m : mask_chans)
ret.img |= m;
}
return ret;
}
Image ConstImage::createSingleChannelMaskFromSet(IntervalSet const& channelSet, bool useAnd) const {
return createSingleChannelMaskFromSet(std::vector<IntervalSet>{channelSet}, useAnd);
}
Image ConstImage::createMultiChannelMaskFromSet(std::vector<IntervalSet> const& channelSets) const { std::vector<IntervalSet> sets = fixBounds(channelSets, type());
std::vector<cv::Mat> mask_chans = getMaskVec(*this, sets);
Image ret;
cv::merge(mask_chans, ret.img);
return ret;
}
Image ConstImage::createMultiChannelMaskFromSet(IntervalSet const& channelSet) const {
return createMultiChannelMaskFromSet(std::vector<IntervalSet>{channelSet});
}
Image ConstImage::convertTo(Type t) const {
Image ret;
img.convertTo(ret.img, toCVType(getBaseType(t)));
return ret;
}
namespace {
template<Type imfu_src_type>
struct ConvertToGrayFunctor {
ConstImage const& src;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), imfu_dst_type};
constexpr bool doRound = isIntegerType(imfu_dst_type);
constexpr double scale = getImageRangeMax(imfu_dst_type) / getImageRangeMax(imfu_src_type);
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
double val = 0.299 * src.at<stype>(x, y, srcChans[0])
+ 0.587 * src.at<stype>(x, y, srcChans[1])
+ 0.114 * src.at<stype>(x, y, srcChans[2]);
dst.at<dtype>(x, y) = static_cast<dtype>((doRound ? 0.5 : 0) + scale * val);
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertToNDIFunctor {
ConstImage const& src;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), imfu_dst_type};
constexpr double scale = getImageRangeMax(imfu_dst_type) / (std::is_signed<dtype>::value ? 1 : 2);
constexpr double offset = std::is_signed<dtype>::value ? 0 : scale;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const stype val1 = src.at<stype>(x, y, srcChans[0]);
const stype val2 = src.at<stype>(x, y, srcChans[1]);
const stype sum = val1 + val2;
const double res = sum == 0 ? 0 : (val1 - val2) * (scale / sum) + offset; // (scale / sum) converts to double
dst.at<dtype>(x, y) = cv::saturate_cast<dtype>(res);
}
}
return dst;
}
};
template<Type imfu_src_type>struct ConvertToBUFunctor {
ConstImage const& src;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), imfu_dst_type};
constexpr double scale = getImageRangeMax(imfu_dst_type) / (std::is_signed<dtype>::value ? 1 : 2);
constexpr double offset = std::is_signed<dtype>::value ? 0 : scale;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const stype red = src.at<stype>(x, y, srcChans[0]);
const stype nir = src.at<stype>(x, y, srcChans[1]);
const stype swir = src.at<stype>(x, y, srcChans[2]);
const double sum_sn = swir + nir;
const double sum_nr = nir + red;
const double res = ( (sum_sn == 0 ? 0 : (swir - nir) / sum_sn)
- (sum_nr == 0 ? 0 : (nir - red) / sum_nr))
* (scale / 2) + offset;
dst.at<dtype>(x, y) = cv::saturate_cast<dtype>(res); // casting non fitting results to integers is UB (wrong number or crash?)
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertToTesseledCapFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double scale = getImageRangeMax(imfu_dst_type) / getImageRangeMax(imfu_src_type) / (std::is_signed<dtype>::value ? 1 : 2);
constexpr double offset = std::is_signed<dtype>::value ? 0 : getImageRangeMax(imfu_dst_type) / 2;
std::array<double, 3*7> f;
if (map == ColorMapping::Landsat_to_TasseledCap) {
constexpr double max = 1 / 2.3103; // if the pixels are all one, this is one over the sum of brightness
f = std::array<double, 3*7>{ // values from: "A Physically-Based Transformation of Thematic Mapper Data - The TM Tasseled Cap" by Crist and Cicone, 1984
// Blue, Green, Red, NIR, SWIR1, SWIR2
+0.3037 * max, +0.2793 * max, +0.4743 * max, +0.5585 * max, +0.5082 * max, +0.1863 * max, 0, // Brightness // TODO: Check if values are correct, for which Landsat: 4, 5, 7, Digital Numbers, Reflectance factors??
-0.2848 * max, -0.2435 * max, -0.5436 * max, +0.7243 * max, +0.0840 * max, -0.1800 * max, 0, // Greenness
+0.1509 * max, +0.1973 * max, +0.3279 * max, +0.3406 * max, -0.7112 * max, -0.4572 * max, 0 // Wetness
};
}
else if (map == ColorMapping::Modis_to_TasseledCap) {
constexpr double max = 1 / 2.6206; // if the pixels are all one, this is one over the sum of brightness
f = std::array<double, 3*7>{ // values from: "MODIS Tasseled Cap Transformation and its Utility" by Zhang et al, 2016
// Red Near-IR Blue Green M-IR M-IR M-IR
// 620-670 841-876 459-479 545-565 1230-1250 1628-1652 2105-2155
+0.3956 * max, +0.4718 * max, +0.3354 * max, +0.3834 * max, +0.3946 * max, +0.3434 * max, +0.2964 * max, // Brightness
-0.3399 * max, +0.5952 * max, -0.2129 * max, -0.2222 * max, +0.4617 * max, -0.1037 * max, -0.4600 * max, // Greenness
+0.10839* max, +0.0912 * max, +0.5065 * max, +0.4040 * max, -0.2410 * max, -0.4658 * max, -0.5306 * max // Wetness
};
}
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
double brightness = 0;
double greeness = 0;
double wetness = 0;
for (unsigned int i = 0; i < srcChans.size(); ++i) {
double src_val = src.at<stype>(x, y, srcChans[i]);
brightness += f[i] * src_val;
greeness += f[i+7] * src_val; wetness += f[i+14] * src_val;
}
brightness = scale * brightness + offset;
greeness = scale * greeness + offset;
wetness = scale * wetness + offset;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(brightness);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(greeness);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(wetness);
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertToYCbCrFunctor {
ConstImage const& src;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double scale = getImageRangeMax(imfu_dst_type) / getImageRangeMax(imfu_src_type);
constexpr double rangeScale = std::is_signed<dtype>::value ? 2 : 1;
constexpr double offset = std::is_signed<dtype>::value ? 0 : getImageRangeMax(imfu_dst_type) / 2;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const double R = src.at<stype>(x, y, srcChans[0]) * scale;
const double G = src.at<stype>(x, y, srcChans[1]) * scale;
const double B = src.at<stype>(x, y, srcChans[2]) * scale;
const double Y = 0.299 * R + 0.587 * G + 0.114 * B;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(Y);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>((B - Y) * (rangeScale / 1.772) + offset);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>((R - Y) * (rangeScale / 1.402) + offset);
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertFromYCbCrFunctor {
ConstImage const& src;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double scale = getImageRangeMax(imfu_dst_type) / getImageRangeMax(imfu_src_type);
constexpr double rangeScale = std::is_signed<stype>::value ? 2 : 1;
constexpr double offset = std::is_signed<stype>::value ? 0 : getImageRangeMax(imfu_src_type) / 2;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const double Y = src.at<stype>(x, y, srcChans[0]) * scale;
const double BY = (src.at<stype>(x, y, srcChans[1]) - offset) * (1.772 * scale / rangeScale);
const double RY = (src.at<stype>(x, y, srcChans[2]) - offset) * (1.402 * scale / rangeScale);
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(Y + RY);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(Y - (0.114 * BY + 0.299 * RY) * (1. / 0.587));
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(Y + BY);
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertToCIEFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double unscale = 1. / getImageRangeMax(imfu_src_type);
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const double R = src.at<stype>(x, y, srcChans[0]) * unscale;
const double G = src.at<stype>(x, y, srcChans[1]) * unscale;
const double B = src.at<stype>(x, y, srcChans[2]) * unscale;
const double X = 0.4339532647955942 * R + 0.3762192150477565 * G + 0.1898275201566493 * B;
const double Y = 0.2126713120163335 * R + 0.7151598465029813 * G + 0.0721688414806854 * B;
const double Z = 0.0177577409665738 * R + 0.1094769889865922 * G + 0.8727652700468339 * B;
constexpr double scale = getImageRangeMax(imfu_dst_type);
if (map == ColorMapping::RGB_to_Lab || map == ColorMapping::RGB_to_Luv) {
constexpr double delta = 6. / 29.;
auto f = [](double t) { return t > std::pow(delta, 3)
? std::cbrt(t)
: t * (1. / (3 * delta * delta)) + 4. / 29.; };
const double L = 1.16 * f(Y) - 0.16; // range [0,1]
double val1;
double val2;
if (map == ColorMapping::RGB_to_Lab) {
constexpr double normScale = !isIntegerType(imfu_dst_type) ? 1 : (std::is_signed<dtype>::value ? 107.862 * scale / (scale + 1) : 206.0972);
constexpr double normOffset = std::is_signed<dtype>::value ? 0 : 0.52335499948568;
val1 /*a*/ = (500 / normScale) * (f(X) - f(Y)) + normOffset; // float / original range [-86.1813, 98.2352], signed int range (scale+1)/scale*[-0.798995939, 0.9107489], unsigned int range [0.10519648, 1]
val2 /*b*/ = (200 / normScale) * (f(Y) - f(Z)) + normOffset; // float / original range [-107.862, 94.4758], signed int range (scale+1)/scale*[-1, 0.8758951], unsigned int range [0, 0.981759]
}
else {
constexpr double normScale = !isIntegerType(imfu_dst_type) ? 1 : (std::is_signed<dtype>::value ? 187.66 : 313.204);
constexpr double normOffset = std::is_signed<dtype>::value ? 0 : 0.400837792620784;
double d = 1. / (X + 15 * Y + 3 * Z + 1e-100);
val1 /*u*/ = (1300 / normScale) * L * (4 * X * d - 0.2009) + normOffset; // float / original range [-78.9986, 187.66], signed int range [-0.4209666, 1], unsigned int range [0.14861049, 1]
val2 /*v*/ = (1300 / normScale) * L * (9 * Y * d - 0.461) + normOffset; // float / original range [-125.544, 116.356], signed int range [-0.668997, 0.620036], unsigned int range [0, 0.77234007]
}
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(L * scale);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(val1 * scale);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(val2 * scale);
}
else { // just to XYZ
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(X * scale);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(Y * scale);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(Z * scale);
}
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertFromCIEFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() { using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
constexpr double scale = getImageRangeMax(imfu_dst_type);
constexpr double unscale = 1. / getImageRangeMax(imfu_src_type);
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
double X = src.at<stype>(x, y, srcChans[0]) * unscale; // might be X or L
double Y = src.at<stype>(x, y, srcChans[1]) * unscale; // might be Y, a or u
double Z = src.at<stype>(x, y, srcChans[2]) * unscale; // might be Z, b or v
if (map == ColorMapping::Lab_to_RGB || map == ColorMapping::Luv_to_RGB) {
constexpr double delta = 6. / 29.;
auto f = [](double t) { return t > delta
? t * t * t
: 3 * delta * delta * (t - 4. / 29.); };
if (map == ColorMapping::Lab_to_RGB) {
const double L = (X + 0.16) * (1. / 1.16);
constexpr double normScale = !isIntegerType(imfu_src_type) ? 1 : (std::is_signed<stype>::value ? 107.862 / unscale / (1/unscale + 1) : 206.0972);
constexpr double normOffset = std::is_signed<stype>::value ? 0 : 0.52335499948568;
const double a = (normScale / 500) * (Y - normOffset);
const double b = (normScale / 200) * (Z - normOffset);
X = f(L + a);
Y = f(L);
Z = f(L - b);
}
else { // map == ColorMapping::Luv_to_RGB
const double& L = X;
if (L == 0) {
Y = Z = 0;
}
else {
constexpr double normScale = !isIntegerType(imfu_src_type) ? 1 : (std::is_signed<stype>::value ? 187.66 : 313.204);
constexpr double normOffset = std::is_signed<stype>::value ? 0 : 0.400837792620784;
double fac = (normScale / 1300) / L;
const double u = fac * (Y - normOffset) + 0.2009;
const double v = fac * (Z - normOffset) + 0.461;
Y = f((L + 0.16) * (1. / 1.16));
X = Y * (9. / 4.) * u / v;
Z = Y * (12 - 3 * u - 20 * v) / (4 * v);
}
}
}
const double R = 3.0799307362950428 * X - 1.5371486218345551 * Y - 0.5427821144604876 * Z;
const double G = -0.9212345908547435 * X + 1.8759903180383508 * Y + 0.0452442728163921 * Z;
const double B = 0.0528909416210833 * X - 0.2040428170607866 * Y + 1.1511518754397034 * Z;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(R * scale);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(G * scale);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(B * scale);
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertToHueFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double unscale = 1. / getImageRangeMax(imfu_src_type);
constexpr double scale = getImageRangeMax(imfu_dst_type); constexpr double Hscale = isIntegerType(imfu_dst_type) ? 1 : 360;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
const double R = src.at<stype>(x, y, srcChans[0]) * unscale;
const double G = src.at<stype>(x, y, srcChans[1]) * unscale;
const double B = src.at<stype>(x, y, srcChans[2]) * unscale;
const double Vmax = std::max({R, G, B});
const double Vmin = std::min({R, G, B});
const double f = (Hscale / 6) / (Vmax - Vmin);
double H;
if (Vmax == Vmin)
H = 0;
else if (Vmax == R)
H = (G - B) * f;
else if (Vmax == G)
H = (B - R) * f + 2 * Hscale / 6;
else //if (Vmax == B)
H = (R - G) * f + 4 * Hscale / 6;
if (H < 0)
H += Hscale;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(scale * H);
if (map == ColorMapping::RGB_to_HSV) {
const double S = Vmax == 0 ? 0 : 1 - Vmin / Vmax;
dst.at<dtype>(x, y, 1) /*S*/ = cv::saturate_cast<dtype>(scale * S);
dst.at<dtype>(x, y, 2) /*V*/ = cv::saturate_cast<dtype>(scale * Vmax);
}
else {
const double L = (Vmin + Vmax) * 0.5;
const double S = L == 0 || L == 1 ? 0 : (Vmax - Vmin) * 0.5 / (L < 0.5 ? L : 1 - L);
dst.at<dtype>(x, y, 1) /*L*/ = cv::saturate_cast<dtype>(scale * L);
dst.at<dtype>(x, y, 2) /*S*/ = cv::saturate_cast<dtype>(scale * S);
}
}
}
return dst;
}
};
template<Type imfu_src_type>
struct ConvertFromHueFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
template<Type imfu_dst_type>
Image operator()() {
using stype = typename DataType<imfu_src_type>::base_type;
using dtype = typename DataType<imfu_dst_type>::base_type;
Image dst{src.size(), getFullType(imfu_dst_type, 3)};
constexpr double scale = getImageRangeMax(imfu_dst_type);
constexpr double unscale = 1. / getImageRangeMax(imfu_src_type);
constexpr double Hscale = isIntegerType(imfu_src_type) ? 1 : 360;
for (int y = 0; y < src.height(); ++y) {
for (int x = 0; x < src.width(); ++x) {
double C;
double m;
if (map == ColorMapping::HSV_to_RGB) {
const double S = src.at<stype>(x, y, srcChans[1]) * unscale;
double V = src.at<stype>(x, y, srcChans[2]) * unscale;
if (S == 0) {
V *= scale;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(V);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(V);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(V);
continue;
}
C = V * S; m = V - C;
}
else {
double L = src.at<stype>(x, y, srcChans[1]) * unscale;
const double S = src.at<stype>(x, y, srcChans[2]) * unscale;
if (S == 0) {
L *= scale;
dst.at<dtype>(x, y, 0) = cv::saturate_cast<dtype>(L);
dst.at<dtype>(x, y, 1) = cv::saturate_cast<dtype>(L);
dst.at<dtype>(x, y, 2) = cv::saturate_cast<dtype>(L);
continue;
}
C = (1 - std::abs(2 * L - 1)) * S;
m = L - C * 0.5;
}
double H = src.at<stype>(x, y, srcChans[0]) * (6 * unscale / Hscale);
while (H < 0) H += 6;
while (H >= 6) H -= 6;
int HCase = static_cast<int>(H);
H -= HCase & (~1); // subtract next lower even number. Same as H = H % 2, but for double, i. e. H = std::fmod(H, 2).
const double X = C * (1 - std::abs(H - 1)) + m;
C += m;
// order for C, X, m
std::array<std::array<int, 3>, 6> order{{{0, 1, 2}, // H in [ 0°, 60°) ==> R = C, G = X, B = m
{1, 0, 2}, // H in [ 60°, 120°) ==> R = X, G = C, B = m
{1, 2, 0}, // H in [120°, 180°) ==> R = m, G = C, B = X
{2, 1, 0}, // H in [180°, 240°) ==> R = m, G = X, B = C
{2, 0, 1}, // H in [240°, 300°) ==> R = X, G = 0, B = C
{0, 2, 1}}};// H in [300°, 360°) ==> R = C, G = 0, B = X
dst.at<dtype>(x, y, order[HCase][0]) = cv::saturate_cast<dtype>(scale * C);
dst.at<dtype>(x, y, order[HCase][1]) = cv::saturate_cast<dtype>(scale * X);
dst.at<dtype>(x, y, order[HCase][2]) = cv::saturate_cast<dtype>(scale * m);
}
}
return dst;
}
};
struct ConvertFunctor {
ConstImage const& src;
ColorMapping map;
std::vector<unsigned int> srcChans;
Type dstType;
template<Type srcType>
Image operator()() {
// using stype = typename DataType<srcType>::base_type;
// Image posMask;
// if (std::numeric_limits<stype>::is_signed) {
// posMask = Image{src.size(), Type::uint8x1};
// posMask.set(0);
// int chans = srcChans.size();
// for (int y = 0; y < src.height(); ++y)
// for (int x = 0; x < src.width(); ++x)
// for (int c = 0; c < chans; ++c)
// posMask.setBoolAt(x, y, 0, posMask.boolAt(x, y, 0) || src.at<stype>(x, y, srcChans[c]) < 0);
// }
switch (map) {
case ColorMapping::RGB_to_Gray:
return CallBaseTypeFunctor::run(ConvertToGrayFunctor<srcType>{src, std::move(srcChans)}, dstType);
case ColorMapping::RGB_to_YCbCr:
return CallBaseTypeFunctor::run(ConvertToYCbCrFunctor<srcType>{src, std::move(srcChans)}, dstType);
case ColorMapping::YCbCr_to_RGB: return CallBaseTypeFunctor::run(ConvertFromYCbCrFunctor<srcType>{src, std::move(srcChans)}, dstType);
case ColorMapping::RGB_to_XYZ:
case ColorMapping::RGB_to_Lab:
case ColorMapping::RGB_to_Luv:
return CallBaseTypeFunctor::run(ConvertToCIEFunctor<srcType>{src, map, std::move(srcChans)}, dstType);
case ColorMapping::XYZ_to_RGB:
case ColorMapping::Lab_to_RGB:
case ColorMapping::Luv_to_RGB:
return CallBaseTypeFunctor::run(ConvertFromCIEFunctor<srcType>{src, map, std::move(srcChans)}, dstType);
case ColorMapping::RGB_to_HSV:
case ColorMapping::RGB_to_HLS:
return CallBaseTypeFunctor::run(ConvertToHueFunctor<srcType>{src, map, std::move(srcChans)}, dstType);
case ColorMapping::HSV_to_RGB:
case ColorMapping::HLS_to_RGB:
return CallBaseTypeFunctor::run(ConvertFromHueFunctor<srcType>{src, map, std::move(srcChans)}, dstType);
case ColorMapping::Pos_Neg_to_NDI:
return CallBaseTypeFunctor::run(ConvertToNDIFunctor<srcType>{src, std::move(srcChans)}, dstType);
case ColorMapping::Red_NIR_SWIR_to_BU:
return CallBaseTypeFunctor::run(ConvertToBUFunctor<srcType>{src, std::move(srcChans)}, dstType);
case ColorMapping::Landsat_to_TasseledCap: // TODO: Add more
case ColorMapping::Modis_to_TasseledCap:
return CallBaseTypeFunctor::run(ConvertToTesseledCapFunctor<srcType>{src, map, std::move(srcChans)}, dstType);
default:
assert(false && "This should not happen. Fix ConvertFunctor!");
return Image{};
}
}
};
constexpr int requiredChannels(ColorMapping map) {
return map == ColorMapping::Pos_Neg_to_NDI ? 2 :
map == ColorMapping::Modis_to_TasseledCap ? 7 :
map == ColorMapping::Landsat_to_TasseledCap ? 6 : // TODO: add more landsat transforms
3;
}
} /* anonymous namespace */
Image ConstImage::convertColor(ColorMapping map, Type result, std::vector<unsigned int> sourceChannels) const {
// using CM = ColorMapping;
for (unsigned int sc : sourceChannels)
if (sc >= channels())
IF_THROW_EXCEPTION(invalid_argument_error("Source channels must be in the range [0, channels-1]. One is: " + std::to_string(sc)));
unsigned int requiredSrcChannels = requiredChannels(map);
if (channels() < requiredSrcChannels)
IF_THROW_EXCEPTION(image_type_error("Color space conversion from " + getFromString(map) + " requires the image to have at least " + std::to_string(requiredSrcChannels) + " channels. The provided image has " + std::to_string(channels()) + " channels."))
<< errinfo_image_type(type());
if (!sourceChannels.empty() && sourceChannels.size() != requiredSrcChannels) {
std::string chans;
for (int sc : sourceChannels)
chans += std::to_string(sc) + ", ";
chans.pop_back();
chans.pop_back();
IF_THROW_EXCEPTION(invalid_argument_error("You provided an invalid number of source channels for a color space conversion from " + getFromString(map)
+ ". Either no source channels should be given or " + std::to_string(requiredSrcChannels)
+ ". You gave: " + std::to_string(sourceChannels.size()) + " [" + chans + "]"));
}
result = getBaseType(result);
if (result == Type::invalid)
result = basetype();
// fill sourceChannel with 0, 1, ...
if (sourceChannels.empty()) {
sourceChannels.resize(requiredSrcChannels);
std::iota(sourceChannels.begin(), sourceChannels.end(), 0);
}
return CallBaseTypeFunctor::run(ConvertFunctor{*this, map, std::move(sourceChannels), result}, basetype());
}
} /* namespace imagefusion */