Read, save and filter images based on GDAL
Currently supports: image filtering of custom operators, image reading, saving, channel separation and other functions
#pragma once #ifdef __GDAL_UTILITY__ #include <gdal/gdal_priv.h> #include <gdal/gdal.h> #include <iostream> #include <string> #include <sstream> #include <fstream> #include <algorithm> #include <list> #include <vector> #include <array> #include <cassert> #include <memory> #include <cstring> #include <initializer_list> #include <random> #include <thread> namespace GD {<!-- --> using namespace std; void init() {<!-- --> GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); } template<typename T> T sum(vector<T> const & amp; v) {<!-- --> T sum = 0; std::for_each(v.begin(), v.end(), [ & amp;](T const & amp; e) {<!-- -->sum + = e; }); return sum; } template<typename T> int partition(vector<T> & amp; vec, int left, int right) {<!-- --> int i = left, j = right + 1; T pivot = vec[left]; while (1) {<!-- --> while (vec[ + + i] < pivot) {<!-- --> if (i == right) break; } while (vec[--j] > pivot) {<!-- --> if (j == left) break; } if (i >= j) break; swap(vec[i], vec[j]); } swap(vec[left], vec[j]); return j; } template<typename T> T qSelect(vector<T> & amp; vec, size_t k) {<!-- --> int low = 0, high = vec.size() - 1; while (low <= high) {<!-- --> int j = partition(vec, low, high); if (j == k) {<!-- --> return vec[k]; } else if (j > k) {<!-- --> high = j - 1; } else if (j < k) {<!-- --> low = j + 1; } } return vec[k]; } template<typename_Ty> struct Matrix {<!-- --> public: static_assert(std::is_arithmetic<_Ty>::value, "_Ty must be numeric type"); //data vector<_Ty*> Data; size_t m_width, m_height; public: Matrix(size_t width = 1, size_t height = 1) :m_width(width), m_height(height) {<!-- --> Data = vector<_Ty*>(height, nullptr); for (int i = 0; i < height; + + i) {<!-- --> Data[i] = new _Ty[width]; memset((void*)Data[i], 0, width * sizeof(_Ty)); } int centerx = width / 2, centery = height / 2; Data[centerx][centery] = 1; } Matrix(std::initializer_list<std::initializer_list<_Ty>> const & amp; init_list) {<!-- --> size_t listW = init_list.begin()->size(), listH = init_list.size(); m_width = listW; m_height = listH; if (listW != m_width || listH != m_height) {<!-- --> cerr << "Matrix init error" << endl; exit(1); } int i = 0; Data = vector<_Ty*>(m_height); for (auto const & amp; list : init_list) {<!-- --> if (list.size() != m_width) {<!-- --> cerr << "Matrix init error" << endl; exit(1); } int j = 0; Data[i] = new _Ty[m_width]; memset((void*)Data[i], 0, m_width * sizeof(_Ty)); for (auto e : list) {<!-- --> Data[i][j + + ] = e; } + + i; } } }; template<typename T> Matrix<T> AverageKernel(size_t size) {<!-- --> if ((size & amp; 1) == 0) {<!-- --> std::cout << "size must be odd" << endl; } Matrix<T> res(size, size); for (int i = 0; i < size; + + i) {<!-- --> for (int j = 0; j < size; j + + ) {<!-- --> res.Data[i][j] = 1.0 / (size * size); } } return res; } template<typename T> Matrix<T> SobelKernelx() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->-1., -2., -1.}, {<!-- --> 0., 0., 0.}, {<!-- --> 1., 2., 1.} }); } template<typename T> Matrix<T> SobelKernely() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->-1., 0., 1.}, {<!-- -->-2., 0., 2.}, {<!-- -->-1., 0., 1.} }); } /** * { {-1., -2., -1.}, { 0., 0., 0.}, { 1., 2., 1.} }); { {-1., 0., 1.}, {-2., 0., 2.}, {-1., 0., 1.} }); */ template<typename T> Matrix<T> LaplaceKernel_8() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->1, 1, 1}, {<!-- -->1, -8, 1}, {<!-- -->1, 1, 1} } ); } template<typename T> Matrix<T> LaplaceKernel_4() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->0, 1, 0}, {<!-- -->1, -4, 1}, {<!-- -->0, 1, 0} } ); } template<typename T> Matrix<T> LaplaceSharpening_4() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->0, 1, 0}, {<!-- -->1, -5, 1}, {<!-- -->0, 1, 0} } ); } template<typename T> Matrix<T> LaplaceSharpening_8() {<!-- --> return Matrix<T>( {<!-- --> {<!-- -->1, 1, 1}, {<!-- -->1, -9, 1}, {<!-- -->1, 1, 1} } ); } template<typename T> Matrix<T> GaussianKernel(size_t size, T stdX, T stdY) {<!-- --> // Get the Gaussian kernel of size*size Matrix<T> ret(size, size); T sum = 0; for (int i = 0; i < size; + + i) {<!-- --> T x = i - size / 2; for (int j = 0; j < size; + + j) {<!-- --> T y = j - size / 2; ret.Data[i][j] = exp(-(x * x / (2 * stdX * stdX) + y * y / (2 * stdY * stdY))); sum + = ret.Data[i][j]; } } for (int i = 0; i < size; + + i) {<!-- --> for (int j = 0; j < size; j + + ) {<!-- --> ret.Data[i][j] /= sum; } } return ret; } class RasterBase {<!-- --> void* pData; public: RasterBase() = default; RasterBase(const string & amp; filepath) {<!-- -->} RasterBase(const char* filepath) {<!-- -->} virtual size_t width() const {<!-- --> return 0; } virtual size_t height() const {<!-- --> return 0; } virtual size_t bands() const {<!-- --> return 0; } virtual void* data() const {<!-- --> return nullptr; } virtual GDALDataType type() const {<!-- --> return GDT_Unknown; } virtual void init(GDALDataset* pDatasetRead, size_t lWidth, size_t lHeight, size_t nBands, GDALDataType datatype) {<!-- -->}; virtual void save(string const & amp;, string const & amp; DriverName = "GTIFF")const {<!-- -->}; virtual vector<array<size_t, 256>> histgram() const {<!-- --> return vector<array<size_t, 256>>(); } virtual RasterBase* equalizeHist(vector<array<size_t, 256>> histgram) {<!-- --> return nullptr; } virtual RasterBase* equalizeHist() {<!-- --> return nullptr; } virtual RasterBase* channel_extraction(size_t channel) {<!-- --> return nullptr; }; virtual RasterBase* filter(Matrix<float> const & amp; kernel) {<!-- --> return nullptr; } virtual RasterBase* filter(Matrix<double> const & amp; kernel) {<!-- --> return nullptr; } virtual RasterBase* MedianFilter() {<!-- --> return nullptr; } virtual RasterBase* addSalt(size_t num) {<!-- --> return nullptr; } virtual RasterBase* Sobel() {<!-- --> return nullptr; } virtual RasterBase* Laplace_4() {<!-- --> return nullptr; } virtual RasterBase* Laplace_8() {<!-- --> return nullptr; } virtual RasterBase* GaussBlur(size_t size, double stdx, double stdy) {<!-- --> return nullptr; } virtual RasterBase* addGaussNoise(size_t size, float stdy) {<!-- --> return nullptr; } virtual RasterBase* Laplace_4_Sharpening() {<!-- --> return nullptr; } virtual RasterBase* Laplace_8_Sharpening() {<!-- --> return nullptr; } virtual RasterBase* Sobel_Sharpening() {<!-- --> return nullptr; } virtual RasterBase* toGray() {<!-- --> return nullptr; } virtual ~RasterBase() {<!-- -->}; }; template<typename_T> _T* Malloc(size_t width, size_t height, size_t bands) {<!-- --> return (_T*)CPLMalloc(width * height * bands * sizeof(_T)); } using RasterPtr = std::shared_ptr<RasterBase>; template<typename T> class Raster : virtual public RasterBase {<!-- --> private: GDALDataset* pDatasetRead; GDALDriver* pDriver; size_t lWidth, lHeight, nBands; T* pData; GDALDataType datatype; template<typename_T> _T & amp; at(size_t x, size_t y) {<!-- --> return pData[x * width * nBands + y]; } template<typename_Ty> vector<_Ty> conv(long long pos, Matrix<_Ty> const & amp; mat) {<!-- --> //Convolve a single pixel, i is the position in the image, and the original image boundary is filled with 0 vector<_Ty> res; size_t total = lHeight * lWidth; long long newpos = pos % total; size_t x = newpos % lWidth, y = newpos / lWidth; long long n = pos / total; size_t matw = mat.m_width, math = mat.m_height; long long matx = matw / 2, maty = math / 2; for (int j = -maty; j <= maty; + + j) {<!-- --> int yj = y + j; for (int i = -matx; i <= matx; + + i) {<!-- --> int xi = x + i; if (xi < 0 || xi >= lHeight || yj < 0 || yj >= lWidth) {<!-- --> res.push_back(0); continue; } long long index = yj * lWidth + xi + n * total; res.push_back(((_Ty)pData[index]) * ((_Ty)mat.Data[maty + j][matx + i])); } } return res; } template<typename_Ty> RasterBase* Laplace(Matrix<_Ty> const & amp; kernel) {<!-- --> RasterPtr pLaplace{<!-- --> filter(kernel) }; RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, 1, datatype) }; size_t Area = lWidth * lHeight; for (size_t i = 0; i < Area; + + i) {<!-- --> double sum{<!-- --> 0.0 }; for (long long channel = 0; channel < nBands; + + channel) {<!-- --> long long index = i + channel * Area; sum + = ((T*)pLaplace->data())[index]; } sum /= 3; ((T*)ret->data())[i] = (T)sum; } return ret; } public: size_t width() const override {<!-- --> return lWidth; } size_t height() const override {<!-- --> return lHeight; } size_t bands() const override {<!-- --> return nBands; } void* data() const override {<!-- --> return (void*)pData; } GDALDataType type() const override {<!-- --> return datatype; } void init(GDALDataset* pDatasetRead, size_t lWidth, size_t lHeight, size_t nBands, GDALDataType datatype)override {<!-- --> this->pDatasetRead = pDatasetRead; this->lWidth = lWidth; this->lHeight = lHeight; this->nBands = nBands; this->datatype = datatype; T* p(Malloc<T>(lWidth, lHeight, nBands)); if (this->pDatasetRead) pDatasetRead->RasterIO(GF_Read, 0, 0, lWidth, lHeight, (void*)p, lWidth, lHeight, datatype, nBands, NULL, 0, 0, 0);//Read image data this->pData = p; } void save(string const & amp; filepath, string const & amp; DriverName = "GTIFF") const override {<!-- --> GDALDataset* pDatasetSave; auto pDriver = GetGDALDriverManager()->GetDriverByName(DriverName.c_str()); pDatasetSave = pDriver->Create(filepath.c_str(), lWidth, lHeight, nBands, datatype, nullptr); assert(pDatasetSave != nullptr); pDatasetSave->RasterIO(GF_Write, 0, 0, lWidth, lHeight, (void*)pData, lWidth, lHeight, datatype, nBands, NULL, 0, 0, 0);//Create image data GDALClose(pDatasetSave); } vector<array<size_t, 256>> histgram() const override {<!-- --> vector<array<size_t, 256>> ret(nBands); for (int i = 0; i < nBands; + + i) {<!-- --> ret[i].fill(0); } long long total = this->lWidth * this->lHeight; for (int channel = 0; channel < nBands; + + channel) {<!-- --> auto p = pData + channel * total; for (int i = 0; i < total; + + i) {<!-- --> + + ret[channel][p[i]]; } } return ret; } RasterBase* equalizeHist(vector<array<size_t, 256>> histgram) override {<!-- --> RasterBase* ret{<!-- --> new Raster<T> }; ret->init(nullptr, lWidth, lHeight, nBands, datatype); size_t total = lHeight * lWidth; memcpy((T*)ret->data(), (T*)this->data(), nBands * total * sizeof(T)); for (int channel = 0; channel < nBands; + + channel) {<!-- --> T* srcbegin = (T*)this->data() + lHeight * lWidth * channel; T* tarbegin = (T*)ret->data() + lHeight * lWidth * channel; for (int i = 1; i < 256; + + i) {<!-- --> histgram[channel][i] + = histgram[channel][i - 1]; } for (int i = 0; i < total; + + i) {<!-- --> tarbegin[i] = static_cast<T>(255.0 * histgram[channel][srcbegin[i]] / total); } } return ret; } RasterBase* equalizeHist()override {<!-- --> return equalizeHist(histgram()); } RasterBase* channel_extraction(size_t channel) {<!-- --> assert(channel < nBands); RasterBase* ret{<!-- --> new Raster<T> }; ret->init(nullptr, lWidth, lHeight, nBands, datatype); size_t total = lHeight * lWidth; T* srcbegin = pData + lHeight * lWidth * channel; T* tarbegin = (T*)ret->data() + lHeight * lWidth * channel; memset(ret->data(), 0, total * sizeof(T)); for (int i = 0; i < total; + + i) {<!-- --> tarbegin[i] = srcbegin[i]; } return ret; } RasterBase* filter(Matrix<float> const & amp; kernel) {<!-- --> if (kernel.m_width % 2 == 0 || kernel.m_height % 2 == 0) {<!-- --> cerr << "kernel size must be odd" << endl; exit(1); } RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; size_t total = lHeight * lWidth * nBands; for (size_t i = 0; i < total; + + i) {<!-- --> auto res = conv(i, kernel); double s = sum(res); if (s < 0) s = fabs(s); if (s > 255) s = 255; ((T*)ret->data())[i] = (T)(s); } return ret; } RasterBase* filter(Matrix<double> const & amp; kernel) {<!-- --> if (kernel.m_width % 2 == 0 || kernel.m_height % 2 == 0) {<!-- --> cerr << "kernel size must be odd" << endl; exit(1); } RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; size_t total = lHeight * lWidth * nBands; for (size_t i = 0; i < total; + + i) {<!-- --> auto res = conv(i, kernel); double s = sum(res); if (s < 0) s = fabs(s); if (s > 255) s = 255; ((T*)ret->data())[i] = (T)(s); } return ret; } RasterBase* MedianFilter()override {<!-- --> RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; size_t total = lHeight * lWidth * nBands; Matrix<int> kernel{<!-- --> {<!-- -->1, 1, 1} ,{<!-- -->1, 1, 1}, {<!-- -- >1, 1, 1} }; for (size_t i = 0; i < total; + + i) {<!-- --> auto res = conv(i, kernel); ((T*)ret->data())[i] = qSelect(res, res.size() / 2); } return ret; } RasterBase* addSalt(size_t num) {<!-- --> size_t limit = lWidth * lHeight; RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; memcpy(ret->data(), this->data(), limit * nBands); for (long long i = 0; i < num; + + i) {<!-- --> int randomIndex = (rand() * rand()) % limit; bool is_white = rand() % 2; for (int channel = 0; channel < nBands; + + channel) {<!-- --> ((T*)ret->data())[randomIndex + channel * limit] = is_white ? 255 : 0; } } return ret; } RasterBase* Sobel() {<!-- --> size_t total = lWidth * lHeight; RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, 1, datatype) }; RasterPtr Sobelx(filter(SobelKernelx<float>())); RasterPtr Sobely(filter(SobelKernely<float>())); for (long long i = 0; i < total; + + i) {<!-- --> double x{<!-- --> .0 }, y{<!-- --> .0 }; for (int channel = 0; channel < nBands; + + channel) {<!-- --> x + = ((T*)Sobelx->data())[i + channel * total]; y + = ((T*)Sobely->data())[i + channel * total]; } double s = sqrt(x * x + y * y); if (s > 255) s = 255; ((T*)ret->data())[i] = (T)s; } return ret; } RasterBase* GaussBlur(size_t size, double stdx, double stdy)override {<!-- --> return filter(GaussianKernel(size, stdx, stdy)); } RasterBase* addGaussNoise(size_t size, float std) {<!-- --> size_t total = lHeight * lWidth * nBands; RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; memcpy(ret->data(), this->data(), total * sizeof(T)); std::default_random_engine e; std::normal_distribution<float> x(0, std); for (long long i = 0; i < size; + + i) {<!-- --> int randomIndex = (rand() * rand()) % total; double s = ((T*)ret->data())[randomIndex] + x(e); if (s < 0) s = 0; if (s > 255) s = 255; ((T*)ret->data())[randomIndex] = (T)(s); } return ret; } RasterBase* Laplace_4()override {<!-- --> return Laplace(LaplaceKernel_4<float>()); } RasterBase* Laplace_8()override {<!-- --> return Laplace(LaplaceKernel_8<float>()); } RasterBase* Laplace_4_Sharpening() {<!-- --> return filter(LaplaceSharpening_4<float>()); } RasterBase* Laplace_8_Sharpening() {<!-- --> return filter(LaplaceSharpening_8<float>()); } RasterBase* Sobel_Sharpening() {<!-- --> size_t total = lWidth * lHeight * nBands; RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, nBands, datatype) }; RasterPtr Sobelx(filter(SobelKernelx<float>())); RasterPtr Sobely(filter(SobelKernely<float>())); for (long long i = 0; i < total; + + i) {<!-- --> T x = ((T*)Sobelx->data())[i]; T y = ((T*)Sobely->data())[i]; double s = sqrt(x * x + y * y) + pData[i]; if (s < 0) s = fabs(s); if (s > 255) s = 255; ((T*)ret->data())[i] = (T)(s); } return ret; } RasterBase* toGray() override {<!-- --> RasterBase* ret{<!-- --> new Raster<T>(lWidth, lHeight, 1, datatype) }; size_t total = lWidth * lHeight; double weights[] = {<!-- --> 0.299, 0.587, 0.114 }; for (size_t i = 0; i < total; + + i) {<!-- --> double sum{<!-- --> 0.0 }; for (long long channel = 0; channel < nBands; + + channel) {<!-- --> long long index = i + channel * total; sum + = weights[channel] * ((T*)pData)[index]; } ((T*)ret->data())[i] = (T)sum; } return ret; } public: Raster() {<!-- --> pData = nullptr; this->lHeight = -1; this->lHeight = -1; pDatasetRead = nullptr; pDriver = nullptr; } Raster(const string & amp; filepath) :Raster(filepath.c_str()) {<!-- -->}; Raster(const char* filepath) {<!-- --> assert(filepath != nullptr); } Raster(size_t width, size_t height, size_t nBands, GDALDataType datatype) {<!-- --> init(nullptr, width, height, nBands, datatype); } ~Raster() {<!-- --> if (pDatasetRead) GDALClose(pDatasetRead);//Close the data set if(pData) delete[] pData; } }; RasterBase* getByType(GDALDataType datatype) {<!-- --> RasterBase* ret = nullptr; switch (datatype) {<!-- --> case GDT_Byte: ret = new Raster<unsigned char>(); break; case GDT_UInt16: ret = new Raster<unsigned short>(); break; case GDT_Int16: ret = new Raster<short>(); break; case GDT_UInt32: ret = new Raster<unsigned int>(); break; case GDT_Int32: ret = new Raster<int>(); break; case GDT_Float32: ret = new Raster<float>(); break; case GDT_Float64: ret = new Raster<double>(); break; } assert(ret != nullptr); return ret; } RasterBase* open(const string & amp; filepath) {<!-- --> GDALDataset* pDatasetRead{<!-- --> nullptr }; size_t lWidth{<!-- --> 0 }, lHeight{<!-- --> 0 }, nBands{<!-- --> 0 }; GDALDataType datatype{<!-- --> GDT_Unknown }; try {<!-- --> pDatasetRead = (GDALDataset*)GDALOpen(filepath.c_str(), GA_ReadOnly); lWidth = pDatasetRead->GetRasterXSize();//Get the image width lHeight = pDatasetRead->GetRasterYSize();//Get the image height nBands = pDatasetRead->GetRasterCount();//Get the number of image bands datatype = pDatasetRead->GetRasterBand(1)->GetRasterDataType();//Get the image format type } catch (std::exception & amp; e) {<!-- --> cerr << e.what() << endl; exit(1); } RasterBase* ptr(getByType(datatype)); ptr->init(pDatasetRead, lWidth, lHeight, nBands, datatype); return ptr; } RasterBase* operator-(RasterBase const & amp; lhs, RasterBase const & amp; rhs) {<!-- --> assert(lhs.width() == rhs.width() & amp; & amp; lhs.height() == rhs.height() & amp; & amp; lhs.bands() == rhs.bands() ); RasterBase* ret{<!-- --> new Raster<unsigned char>(lhs.width(), lhs.height(), lhs.bands(), lhs.type()) }; size_t total = lhs.width() * lhs.height() * lhs.bands(); for (size_t i = 0; i < total; + + i) {<!-- --> double diff = ((unsigned char*)lhs.data())[i] - ((unsigned char*)rhs.data())[i]; if (diff < 0) diff = 0; ((unsigned char*)ret->data())[i] = (unsigned char)diff; } return ret; } } #endif // __GDAL__UTILITY__H__