From 38e4ef464c16d42c92db268f75db93fc1aabb72b Mon Sep 17 00:00:00 2001 From: Matthew Date: Mon, 30 Dec 2024 23:30:20 +0800 Subject: [PATCH] format --- bin/hdrplus.cpp | 22 +- include/hdrplus/finish.h | 471 +++++----- include/hdrplus/params.h | 98 +-- src/align.cpp | 1764 +++++++++++++++++++------------------- src/bayer_image.cpp | 370 ++++---- src/finish.cpp | 1449 ++++++++++++++++--------------- src/hdrplus_pipeline.cpp | 156 ++-- src/merge.cpp | 636 +++++++------- 8 files changed, 2496 insertions(+), 2470 deletions(-) diff --git a/bin/hdrplus.cpp b/bin/hdrplus.cpp index 948785c..8e4d7a2 100644 --- a/bin/hdrplus.cpp +++ b/bin/hdrplus.cpp @@ -1,19 +1,19 @@ #include "hdrplus/hdrplus_pipeline.h" -int main( int argc, char** argv ) +int main(int argc, char** argv) { - std::vector paths; - for (int idx = 2; idx < argc; idx++) - { - paths.push_back(argv[idx]); - } + std::vector paths; + for (int idx = 2; idx < argc; idx++) + { + paths.push_back(argv[idx]); + } - cv::Mat rgb; + cv::Mat rgb; - hdrplus::hdrplus_pipeline pipeline; - pipeline.run_pipeline( paths, 0, rgb); + hdrplus::hdrplus_pipeline pipeline; + pipeline.run_pipeline(paths, 0, rgb); - cv::imwrite(argv[1], rgb); + cv::imwrite(argv[1], rgb); - return 0; + return 0; } \ No newline at end of file diff --git a/include/hdrplus/finish.h b/include/hdrplus/finish.h index fdc86eb..54dc317 100644 --- a/include/hdrplus/finish.h +++ b/include/hdrplus/finish.h @@ -13,239 +13,242 @@ namespace hdrplus { - uint16_t uGammaCompress_1pix(float x, float threshold,float gainMin,float gainMax,float exponent); - uint16_t uGammaDecompress_1pix(float x, float threshold,float gainMin,float gainMax,float exponent); - cv::Mat uGammaCompress_(cv::Mat m,float threshold,float gainMin,float gainMax,float exponent); - cv::Mat uGammaDecompress_(cv::Mat m,float threshold,float gainMin,float gainMax,float exponent); - cv::Mat gammasRGB(cv::Mat img, bool mode); - - -class finish -{ - public: - cv::Mat mergedBayer; // merged image from Merge Module - std::string burstPath; // path to burst images - std::vector rawPathList; // a list or array of the path to all burst imgs under burst Path - int refIdx; // index of the reference img - Parameters params; - cv::Mat rawReference; - // LibRaw libraw_processor_finish; - bayer_image* refBayer; - - std::string mergedImgPath; - finish() - { - refBayer = NULL; - } - -// please use this initialization after merging part finish - finish(std::string burstPath, cv::Mat mergedBayer,int refIdx) { - refBayer = NULL; - this->refIdx = refIdx; - this->burstPath = burstPath; - this->mergedBayer = mergedBayer; - } - -// for local testing only - finish(std::string burstPath, std::string mergedBayerPath,int refIdx){ - this->refIdx = refIdx; - this->burstPath = burstPath; - this->mergedBayer = loadFromCSV(mergedBayerPath, CV_16UC1);// - load_rawPathList(burstPath); - refBayer= new bayer_image(this->rawPathList[refIdx]); - this->rawReference = refBayer->raw_image;//;grayscale_image - - // initialize parameters in libraw_processor_finish - setLibRawParams(); - showParams(); - - std::cout<<"Finish init() finished!"<& libraw_ptr, cv::Mat B); - - // postprocess - // cv::Mat postprocess(std::shared_ptr& libraw_ptr); - - void showImg(cv::Mat img) - { - int ch = CV_MAT_CN(CV_8UC1); - - // cv::Mat tmp(4208,3120,CV_16UC1); - cv::Mat tmp(img); - // uint16_t* ptr_tmp = (uint16_t*)tmp.data; - // uint16_t* ptr_img = (uint16_t*)img.data; - // // col major to row major - // for(int r = 0; r < tmp.rows; r++) { - // for(int c = 0; c < tmp.cols; c++) { - // *(ptr_tmp+r*tmp.cols+c) = *(ptr_img+c*tmp.rows+r)/2048.0*255.0; - // } - // } - // std::cout<<"height="<mergedBayer.size()<rawPathList){ - std::cout< dvals; - std::stringstream ss(line); - std::string val; - // int count=0; - while (getline(ss, val, ',')) - { - dvals.push_back(stod(val));//*255.0/2048.0 - // count++; - } - // std::cout<d_name; // current filepath that ptr points to - if (ptr->d_type != 8 && ptr->d_type != 4) { // not normal file or dir - return; - } - // only need normal files - if (ptr->d_type == 8) { - if (strcmp(ptr->d_name, ".") != 0 && strcmp(ptr->d_name, "..") != 0) { - if (strstr(ptr->d_name, ".dng")) { - rawPathList.emplace_back(sub_file); - } - } - } - } - // close root dir - closedir(pDir); - } - - void setLibRawParams(){ - refBayer->libraw_processor->imgdata.params.user_qual = params.rawpyArgs.demosaic_algorithm; - refBayer->libraw_processor->imgdata.params.half_size = params.rawpyArgs.half_size; - refBayer->libraw_processor->imgdata.params.use_camera_wb = params.rawpyArgs.use_camera_wb; - refBayer->libraw_processor->imgdata.params.use_auto_wb = params.rawpyArgs.use_auto_wb; - refBayer->libraw_processor->imgdata.params.no_auto_bright = params.rawpyArgs.no_auto_bright; - refBayer->libraw_processor->imgdata.params.output_color = params.rawpyArgs.output_color; - refBayer->libraw_processor->imgdata.params.gamm[0] = params.rawpyArgs.gamma[0]; - refBayer->libraw_processor->imgdata.params.gamm[1] = params.rawpyArgs.gamma[1]; - refBayer->libraw_processor->imgdata.params.output_bps = params.rawpyArgs.output_bps; - - // libraw_processor_finish.imgdata.params.user_qual = params.rawpyArgs.demosaic_algorithm; - // libraw_processor_finish.imgdata.params.half_size = params.rawpyArgs.half_size; - // libraw_processor_finish.imgdata.params.use_camera_wb = params.rawpyArgs.use_camera_wb; - // libraw_processor_finish.imgdata.params.use_auto_wb = params.rawpyArgs.use_auto_wb; - // libraw_processor_finish.imgdata.params.no_auto_bright = params.rawpyArgs.no_auto_bright; - // libraw_processor_finish.imgdata.params.output_color = params.rawpyArgs.output_color; - // libraw_processor_finish.imgdata.params.gamm[0] = params.rawpyArgs.gamma[0]; - // libraw_processor_finish.imgdata.params.gamm[1] = params.rawpyArgs.gamma[1]; - // libraw_processor_finish.imgdata.params.output_bps = params.rawpyArgs.output_bps; - } - - -}; + uint16_t uGammaCompress_1pix(float x, float threshold, float gainMin, float gainMax, float exponent); + uint16_t uGammaDecompress_1pix(float x, float threshold, float gainMin, float gainMax, float exponent); + cv::Mat uGammaCompress_(cv::Mat m, float threshold, float gainMin, float gainMax, float exponent); + cv::Mat uGammaDecompress_(cv::Mat m, float threshold, float gainMin, float gainMax, float exponent); + cv::Mat gammasRGB(cv::Mat img, bool mode); + + + class finish + { + public: + cv::Mat mergedBayer; // merged image from Merge Module + std::string burstPath; // path to burst images + std::vector rawPathList; // a list or array of the path to all burst imgs under burst Path + int refIdx; // index of the reference img + Parameters params; + cv::Mat rawReference; + // LibRaw libraw_processor_finish; + bayer_image* refBayer; + + std::string mergedImgPath; + finish() + { + refBayer = NULL; + } + + // please use this initialization after merging part finish + finish(std::string burstPath, cv::Mat mergedBayer, int refIdx) { + refBayer = NULL; + this->refIdx = refIdx; + this->burstPath = burstPath; + this->mergedBayer = mergedBayer; + } + + // for local testing only + finish(std::string burstPath, std::string mergedBayerPath, int refIdx) { + this->refIdx = refIdx; + this->burstPath = burstPath; + this->mergedBayer = loadFromCSV(mergedBayerPath, CV_16UC1);// + load_rawPathList(burstPath); + refBayer = new bayer_image(this->rawPathList[refIdx]); + this->rawReference = refBayer->raw_image;//;grayscale_image + + // initialize parameters in libraw_processor_finish + setLibRawParams(); + showParams(); + + std::cout << "Finish init() finished!" << std::endl; + } + + + + ~finish() + { + if (refBayer != NULL) + { + delete refBayer; + refBayer = NULL; + } + } + + // finish pipeline func + // void process(std::string burstPath, cv::Mat mergedBayer,int refIdx); + void process(const hdrplus::burst& burst_images, cv::Mat& finalOutputImage); + + // replace Mat a with Mat b + void copy_mat_16U(cv::Mat& A, cv::Mat B); + void copy_rawImg2libraw(std::shared_ptr& libraw_ptr, cv::Mat B); + + // postprocess + // cv::Mat postprocess(std::shared_ptr& libraw_ptr); + + void showImg(cv::Mat img) + { + int ch = CV_MAT_CN(CV_8UC1); + + // cv::Mat tmp(4208,3120,CV_16UC1); + cv::Mat tmp(img); + // uint16_t* ptr_tmp = (uint16_t*)tmp.data; + // uint16_t* ptr_img = (uint16_t*)img.data; + // // col major to row major + // for(int r = 0; r < tmp.rows; r++) { + // for(int c = 0; c < tmp.cols; c++) { + // *(ptr_tmp+r*tmp.cols+c) = *(ptr_img+c*tmp.rows+r)/2048.0*255.0; + // } + // } + // std::cout<<"height="<mergedBayer.size() << std::endl; + } + + void showMat(cv::Mat img) { + std::cout << "size=" << img.size() << std::endl; + std::cout << "type=" << img.type() << std::endl; + } + + void showParams() + { + std::cout << "Parameters:" << std::endl; + std::cout << "tuning_ltmGain = " << this->params.tuning.ltmGain << std::endl; + std::cout << "tuning_gtmContrast = " << this->params.tuning.gtmContrast << std::endl; + for (auto key_val : this->params.flags) { + std::cout << key_val.first << "," << key_val.second << std::endl; + } + std::cout << "demosaic_algorithm = " << refBayer->libraw_processor->imgdata.params.user_qual << std::endl; + std::cout << "half_size = " << refBayer->libraw_processor->imgdata.params.half_size << std::endl; + std::cout << "use_camera_wb = " << refBayer->libraw_processor->imgdata.params.use_camera_wb << std::endl; + std::cout << "use_auto_wb = " << refBayer->libraw_processor->imgdata.params.use_auto_wb << std::endl; + std::cout << "no_auto_bright = " << refBayer->libraw_processor->imgdata.params.no_auto_bright << std::endl; + std::cout << "output_color = " << refBayer->libraw_processor->imgdata.params.output_color << std::endl; + std::cout << "gamma[0] = " << refBayer->libraw_processor->imgdata.params.gamm[0] << std::endl; + std::cout << "gamma[1] = " << refBayer->libraw_processor->imgdata.params.gamm[1] << std::endl; + std::cout << "output_bps = " << refBayer->libraw_processor->imgdata.params.output_bps << std::endl; + + // std::cout<<"demosaic_algorithm = "<rawPathList) { + std::cout << pth << std::endl; + } + std::cout << "====================" << std::endl; + } + + + + + private: + cv::Mat loadFromCSV(const std::string& path, int opencv_type) + { + cv::Mat m; + std::ifstream csvFile(path); + + std::string line; + + while (getline(csvFile, line)) + { + std::vector dvals; + std::stringstream ss(line); + std::string val; + // int count=0; + while (getline(ss, val, ',')) + { + dvals.push_back(stod(val));//*255.0/2048.0 + // count++; + } + // std::cout<d_name; // current filepath that ptr points to + if (ptr->d_type != 8 && ptr->d_type != 4) { // not normal file or dir + return; + } + // only need normal files + if (ptr->d_type == 8) { + if (strcmp(ptr->d_name, ".") != 0 && strcmp(ptr->d_name, "..") != 0) { + if (strstr(ptr->d_name, ".dng")) { + rawPathList.emplace_back(sub_file); + } + } + } + } + // close root dir + closedir(pDir); + } + + void setLibRawParams() { + refBayer->libraw_processor->imgdata.params.user_qual = params.rawpyArgs.demosaic_algorithm; + refBayer->libraw_processor->imgdata.params.half_size = params.rawpyArgs.half_size; + refBayer->libraw_processor->imgdata.params.use_camera_wb = params.rawpyArgs.use_camera_wb; + refBayer->libraw_processor->imgdata.params.use_auto_wb = params.rawpyArgs.use_auto_wb; + refBayer->libraw_processor->imgdata.params.no_auto_bright = params.rawpyArgs.no_auto_bright; + refBayer->libraw_processor->imgdata.params.output_color = params.rawpyArgs.output_color; + refBayer->libraw_processor->imgdata.params.gamm[0] = params.rawpyArgs.gamma[0]; + refBayer->libraw_processor->imgdata.params.gamm[1] = params.rawpyArgs.gamma[1]; + refBayer->libraw_processor->imgdata.params.output_bps = params.rawpyArgs.output_bps; + + // libraw_processor_finish.imgdata.params.user_qual = params.rawpyArgs.demosaic_algorithm; + // libraw_processor_finish.imgdata.params.half_size = params.rawpyArgs.half_size; + // libraw_processor_finish.imgdata.params.use_camera_wb = params.rawpyArgs.use_camera_wb; + // libraw_processor_finish.imgdata.params.use_auto_wb = params.rawpyArgs.use_auto_wb; + // libraw_processor_finish.imgdata.params.no_auto_bright = params.rawpyArgs.no_auto_bright; + // libraw_processor_finish.imgdata.params.output_color = params.rawpyArgs.output_color; + // libraw_processor_finish.imgdata.params.gamm[0] = params.rawpyArgs.gamma[0]; + // libraw_processor_finish.imgdata.params.gamm[1] = params.rawpyArgs.gamma[1]; + // libraw_processor_finish.imgdata.params.output_bps = params.rawpyArgs.output_bps; + } + + + }; } // namespace hdrplus diff --git a/include/hdrplus/params.h b/include/hdrplus/params.h index 6bf1744..7408b4d 100644 --- a/include/hdrplus/params.h +++ b/include/hdrplus/params.h @@ -8,62 +8,62 @@ namespace hdrplus { -class RawpyArgs{ - public: - int demosaic_algorithm = 3;// 3 - AHD interpolation <->int user_qual - bool half_size = false; - bool use_camera_wb = true; - bool use_auto_wb = false; - bool no_auto_bright = true; - int output_color = LIBRAW_COLORSPACE_sRGB; - int gamma[2] = {1,1}; //# gamma correction not applied by rawpy (not quite understand) - int output_bps = 16; -}; + class RawpyArgs { + public: + int demosaic_algorithm = 3;// 3 - AHD interpolation <->int user_qual + bool half_size = false; + bool use_camera_wb = true; + bool use_auto_wb = false; + bool no_auto_bright = true; + int output_color = LIBRAW_COLORSPACE_sRGB; + int gamma[2] = { 1,1 }; //# gamma correction not applied by rawpy (not quite understand) + int output_bps = 16; + }; -class Options{ - public: - std::string input = ""; - std::string output = ""; - std::string mode = "full"; //'full' 'align' 'merge' 'finish' - int reference = 0; - float temporalfactor=75.0; - float spatialfactor = 0.1; - int ltmGain=-1; - double gtmContrast=0.075; - int verbose=2; // (0, 1, 2, 3, 4, 5) + class Options { + public: + std::string input = ""; + std::string output = ""; + std::string mode = "full"; //'full' 'align' 'merge' 'finish' + int reference = 0; + float temporalfactor = 75.0; + float spatialfactor = 0.1; + int ltmGain = -1; + double gtmContrast = 0.075; + int verbose = 2; // (0, 1, 2, 3, 4, 5) -}; + }; -class Tuning{ - public: - std::string ltmGain = "auto"; - double gtmContrast = 0.075; - std::vector sharpenAmount{1,0.5,0.5}; - std::vector sharpenSigma{1,2,4}; - std::vector sharpenThreshold{0.02,0.04,0.06}; -}; + class Tuning { + public: + std::string ltmGain = "auto"; + double gtmContrast = 0.075; + std::vector sharpenAmount{ 1,0.5,0.5 }; + std::vector sharpenSigma{ 1,2,4 }; + std::vector sharpenThreshold{ 0.02,0.04,0.06 }; + }; -class Parameters{ - public: - std::unordered_map flags; + class Parameters { + public: + std::unordered_map flags; - RawpyArgs rawpyArgs; - Options options; - Tuning tuning; + RawpyArgs rawpyArgs; + Options options; + Tuning tuning; - Parameters() - { - const char* keys[] = {"writeReferenceImage", "writeGammaReference", "writeMergedImage", "writeGammaMerged", - "writeShortExposure", "writeLongExposure", "writeFusedExposure", "writeLTMImage", - "writeLTMGamma", "writeGTMImage", "writeReferenceFinal", "writeFinalImage"}; - for (int idx = 0; idx < sizeof(keys) / sizeof(const char*); idx++) { - flags[keys[idx]] = true; - } - } + Parameters() + { + const char* keys[] = { "writeReferenceImage", "writeGammaReference", "writeMergedImage", "writeGammaMerged", + "writeShortExposure", "writeLongExposure", "writeFusedExposure", "writeLTMImage", + "writeLTMGamma", "writeGTMImage", "writeReferenceFinal", "writeFinalImage" }; + for (int idx = 0; idx < sizeof(keys) / sizeof(const char*); idx++) { + flags[keys[idx]] = true; + } + } -}; + }; -cv::Mat postprocess(std::shared_ptr& libraw_ptr, RawpyArgs rawpyArgs); -void setParams(std::shared_ptr& libraw_ptr, RawpyArgs rawpyArgs); + cv::Mat postprocess(std::shared_ptr& libraw_ptr, RawpyArgs rawpyArgs); + void setParams(std::shared_ptr& libraw_ptr, RawpyArgs rawpyArgs); } // namespace hdrplus diff --git a/src/align.cpp b/src/align.cpp index 46efa04..390c6cc 100644 --- a/src/align.cpp +++ b/src/align.cpp @@ -13,986 +13,986 @@ namespace hdrplus { -// Function declration -static void build_per_grayimg_pyramid( \ - std::vector& images_pyramid, \ - const cv::Mat& src_image, \ - const std::vector& inv_scale_factors ); - - -template< int pyramid_scale_factor_prev_curr, int tilesize_scale_factor_prev_curr, int tile_size > -static void build_upsampled_prev_aligement( \ - const std::vector>>& src_alignment, \ - std::vector>>& dst_alignment, \ - int num_tiles_h, int num_tiles_w, \ - const cv::Mat& ref_img, const cv::Mat& alt_img, \ - bool consider_nbr = false ); - - -template< typename data_type, typename return_type, int tile_size > -static unsigned long long l1_distance( const cv::Mat& img1, const cv::Mat& img2, \ - int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ - int img2_tile_row_start_idx, int img2_tile_col_start_idx ); + // Function declration + static void build_per_grayimg_pyramid(\ + std::vector& images_pyramid, \ + const cv::Mat& src_image, \ + const std::vector& inv_scale_factors); + + + template< int pyramid_scale_factor_prev_curr, int tilesize_scale_factor_prev_curr, int tile_size > + static void build_upsampled_prev_aligement(\ + const std::vector>>& src_alignment, \ + std::vector>>& dst_alignment, \ + int num_tiles_h, int num_tiles_w, \ + const cv::Mat& ref_img, const cv::Mat& alt_img, \ + bool consider_nbr = false); + + + template< typename data_type, typename return_type, int tile_size > + static unsigned long long l1_distance(const cv::Mat& img1, const cv::Mat& img2, \ + int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ + int img2_tile_row_start_idx, int img2_tile_col_start_idx); + + + template< typename data_type, typename return_type, int tile_size > + static return_type l2_distance(const cv::Mat& img1, const cv::Mat& img2, \ + int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ + int img2_tile_row_start_idx, int img2_tile_col_start_idx); + + + static void align_image_level(\ + const cv::Mat& ref_img, \ + const cv::Mat& alt_img, \ + std::vector>>& prev_aligement, \ + std::vector>>& curr_alignment, \ + int scale_factor_prev_curr, \ + int curr_tile_size, \ + int prev_tile_size, \ + int search_radiou, \ + int distance_type); + + + // Function Implementations + + + // static function only visible within file + static void build_per_grayimg_pyramid(\ + std::vector& images_pyramid, \ + const cv::Mat& src_image, \ + const std::vector& inv_scale_factors) + { + // #ifndef NDEBUG + // printf("%s::%s build_per_grayimg_pyramid start with scale factor : ", __FILE__, __func__ ); + // for ( int i = 0; i < inv_scale_factors.size(); ++i ) + // { + // printf("%d ", inv_scale_factors.at( i )); + // } + // printf("\n"); + // #endif + + images_pyramid.resize(inv_scale_factors.size()); + + for (size_t i = 0; i < inv_scale_factors.size(); ++i) + { + cv::Mat blur_image; + cv::Mat downsample_image; + + switch (inv_scale_factors[i]) + { + case 1: + images_pyramid[i] = src_image.clone(); + // cv::Mat use reference count, will not create deep copy + downsample_image = src_image; + break; + case 2: + // printf("(2) downsample with gaussian sigma %.2f", inv_scale_factors[ i ] * 0.5 ); + // // Gaussian blur + cv::GaussianBlur(images_pyramid.at(i - 1), blur_image, cv::Size(0, 0), inv_scale_factors[i] * 0.5); + + // // Downsample + downsample_image = downsample_nearest_neighbour(blur_image); + // downsample_image = downsample_nearest_neighbour( images_pyramid.at( i-1 ) ); + + // Add + images_pyramid.at(i) = downsample_image.clone(); + + break; + case 4: + // printf("(4) downsample with gaussian sigma %.2f", inv_scale_factors[ i ] * 0.5 ); + cv::GaussianBlur(images_pyramid.at(i - 1), blur_image, cv::Size(0, 0), inv_scale_factors[i] * 0.5); + downsample_image = downsample_nearest_neighbour(blur_image); + // downsample_image = downsample_nearest_neighbour( images_pyramid.at( i-1 ) ); + images_pyramid.at(i) = downsample_image.clone(); + break; + default: +#ifdef __ANDROID__ + break; +#else + throw std::runtime_error("inv scale factor " + std::to_string(inv_scale_factors[i]) + "invalid"); +#endif + } + } + } -template< typename data_type, typename return_type, int tile_size > -static return_type l2_distance( const cv::Mat& img1, const cv::Mat& img2, \ - int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ - int img2_tile_row_start_idx, int img2_tile_col_start_idx ); + static bool operator!=(const std::pair& lhs, const std::pair& rhs) + { + return lhs.first != rhs.first || lhs.second != rhs.second; + } -static void align_image_level( \ - const cv::Mat& ref_img, \ - const cv::Mat& alt_img, \ - std::vector>>& prev_aligement, \ - std::vector>>& curr_alignment, \ - int scale_factor_prev_curr, \ - int curr_tile_size, \ - int prev_tile_size, \ - int search_radiou, \ - int distance_type ); + template< int pyramid_scale_factor_prev_curr, int tilesize_scale_factor_prev_curr, int tile_size > + static void build_upsampled_prev_aligement(\ + const std::vector>>& src_alignment, \ + std::vector>>& dst_alignment, \ + int num_tiles_h, int num_tiles_w, \ + const cv::Mat& ref_img, const cv::Mat& alt_img, \ + bool consider_nbr) + { + int src_num_tiles_h = src_alignment.size(); + int src_num_tiles_w = src_alignment[0].size(); + constexpr int repeat_factor = pyramid_scale_factor_prev_curr / tilesize_scale_factor_prev_curr; -// Function Implementations + // printf("build_upsampled_prev_aligement with scale factor %d, repeat factor %d, tile size factor %d\n", \ + // pyramid_scale_factor_prev_curr, repeat_factor, tilesize_scale_factor_prev_curr ); + int dst_num_tiles_main_h = src_num_tiles_h * repeat_factor; + int dst_num_tiles_main_w = src_num_tiles_w * repeat_factor; -// static function only visible within file -static void build_per_grayimg_pyramid( \ - std::vector& images_pyramid, \ - const cv::Mat& src_image, \ - const std::vector& inv_scale_factors ) -{ - // #ifndef NDEBUG - // printf("%s::%s build_per_grayimg_pyramid start with scale factor : ", __FILE__, __func__ ); - // for ( int i = 0; i < inv_scale_factors.size(); ++i ) - // { - // printf("%d ", inv_scale_factors.at( i )); - // } - // printf("\n"); - // #endif - - images_pyramid.resize( inv_scale_factors.size() ); - - for ( size_t i = 0; i < inv_scale_factors.size(); ++i ) - { - cv::Mat blur_image; - cv::Mat downsample_image; - - switch ( inv_scale_factors[ i ] ) - { - case 1: - images_pyramid[ i ] = src_image.clone(); - // cv::Mat use reference count, will not create deep copy - downsample_image = src_image; - break; - case 2: - // printf("(2) downsample with gaussian sigma %.2f", inv_scale_factors[ i ] * 0.5 ); - // // Gaussian blur - cv::GaussianBlur( images_pyramid.at( i-1 ), blur_image, cv::Size(0, 0), inv_scale_factors[ i ] * 0.5 ); - - // // Downsample - downsample_image = downsample_nearest_neighbour( blur_image ); - // downsample_image = downsample_nearest_neighbour( images_pyramid.at( i-1 ) ); - - // Add - images_pyramid.at( i ) = downsample_image.clone(); - - break; - case 4: - // printf("(4) downsample with gaussian sigma %.2f", inv_scale_factors[ i ] * 0.5 ); - cv::GaussianBlur( images_pyramid.at( i-1 ), blur_image, cv::Size(0, 0), inv_scale_factors[ i ] * 0.5 ); - downsample_image = downsample_nearest_neighbour( blur_image ); - // downsample_image = downsample_nearest_neighbour( images_pyramid.at( i-1 ) ); - images_pyramid.at( i ) = downsample_image.clone(); - break; - default: + if (dst_num_tiles_main_h > num_tiles_h || dst_num_tiles_main_w > num_tiles_w) + { #ifdef __ANDROID__ - break; + return; #else - throw std::runtime_error("inv scale factor " + std::to_string( inv_scale_factors[ i ]) + "invalid" ); + throw std::runtime_error("current level number of tiles smaller than upsampled tiles\n"); +#endif + } + + // Allocate data for dst_alignment + // NOTE: number of tiles h, number of tiles w might be different from dst_num_tiles_main_h, dst_num_tiles_main_w + // For tiles between num_tile_h and dst_num_tiles_main_h, use (0,0) + dst_alignment.resize(num_tiles_h, std::vector>(num_tiles_w, std::pair(0, 0))); + + // Upsample alignment +#pragma omp parallel for collapse(2) + for (int row_i = 0; row_i < src_num_tiles_h; row_i++) + { + for (int col_i = 0; col_i < src_num_tiles_w; col_i++) + { + // Scale alignment + std::pair align_i = src_alignment[row_i][col_i]; + align_i.first *= pyramid_scale_factor_prev_curr; + align_i.second *= pyramid_scale_factor_prev_curr; + + // repeat + UNROLL_LOOP(repeat_factor) + for (int repeat_row_i = 0; repeat_row_i < repeat_factor; ++repeat_row_i) + { + int repeat_row_i_offset = row_i * repeat_factor + repeat_row_i; + UNROLL_LOOP(repeat_factor) + for (int repeat_col_i = 0; repeat_col_i < repeat_factor; ++repeat_col_i) + { + int repeat_col_i_offset = col_i * repeat_factor + repeat_col_i; + dst_alignment[repeat_row_i_offset][repeat_col_i_offset] = align_i; + } + } + } + } + + if (consider_nbr) + { + // Copy consurtctor + std::vector>> upsampled_alignment{ dst_alignment }; + + // Distance function + unsigned long long(*distance_func_ptr)(const cv::Mat&, const cv::Mat&, int, int, int, int) = \ + &l1_distance; + +#pragma omp parallel for collapse(2) + for (int tile_row_i = 0; tile_row_i < num_tiles_h; tile_row_i++) + { + for (int tile_col_i = 0; tile_col_i < num_tiles_w; tile_col_i++) + { + const auto& curr_align_i = upsampled_alignment[tile_row_i][tile_col_i]; + + // Container for nbr alignment pair + std::vector> nbrs_align_i; + + // Consider 4 neighbour's alignment + // Only compute distance if alignment is different + if (tile_col_i > 0) + { + const auto& nbr1_align_i = upsampled_alignment[tile_row_i + 0][tile_col_i - 1]; + if (curr_align_i != nbr1_align_i) nbrs_align_i.emplace_back(nbr1_align_i); + } + + if (tile_col_i < num_tiles_w - 1) + { + const auto& nbr2_align_i = upsampled_alignment[tile_row_i + 0][tile_col_i + 1]; + if (curr_align_i != nbr2_align_i) nbrs_align_i.emplace_back(nbr2_align_i); + } + + if (tile_row_i > 0) + { + const auto& nbr3_align_i = upsampled_alignment[tile_row_i - 1][tile_col_i + 0]; + if (curr_align_i != nbr3_align_i) nbrs_align_i.emplace_back(nbr3_align_i); + } + + if (tile_row_i < num_tiles_h - 1) + { + const auto& nbr4_align_i = upsampled_alignment[tile_row_i + 1][tile_col_i + 0]; + if (curr_align_i != nbr4_align_i) nbrs_align_i.emplace_back(nbr4_align_i); + } + + // If there is a nbr alignment that need to be considered. Compute distance + if (!nbrs_align_i.empty()) + { + int ref_tile_row_start_idx_i = tile_row_i * tile_size / 2; + int ref_tile_col_start_idx_i = tile_col_i * tile_size / 2; + + // curr_align_i's distance + auto curr_align_i_distance = distance_func_ptr( + ref_img, alt_img, \ + ref_tile_row_start_idx_i, \ + ref_tile_col_start_idx_i, \ + ref_tile_row_start_idx_i + curr_align_i.first, \ + ref_tile_col_start_idx_i + curr_align_i.second); + + for (const auto& nbr_align_i : nbrs_align_i) + { + auto nbr_align_i_distance = distance_func_ptr( + ref_img, alt_img, \ + ref_tile_row_start_idx_i, \ + ref_tile_col_start_idx_i, \ + ref_tile_row_start_idx_i + nbr_align_i.first, \ + ref_tile_col_start_idx_i + nbr_align_i.second); + + if (nbr_align_i_distance < curr_align_i_distance) + { +#ifdef NDEBUG + printf("tile [%d, %d] update align, prev align (%d, %d) curr align (%d, %d), prev distance %d curr distance %d\n", \ + tile_row_i, tile_col_i, \ + curr_align_i.first, curr_align_i.second, \ + nbr_align_i.first, nbr_align_i.second, \ + int(curr_align_i_distance), int(nbr_align_i_distance)); #endif - } - } -} - - -static bool operator!=( const std::pair& lhs, const std::pair& rhs ) -{ - return lhs.first != rhs.first || lhs.second != rhs.second; -} - - -template< int pyramid_scale_factor_prev_curr, int tilesize_scale_factor_prev_curr, int tile_size > -static void build_upsampled_prev_aligement( \ - const std::vector>>& src_alignment, \ - std::vector>>& dst_alignment, \ - int num_tiles_h, int num_tiles_w, \ - const cv::Mat& ref_img, const cv::Mat& alt_img, \ - bool consider_nbr ) -{ - int src_num_tiles_h = src_alignment.size(); - int src_num_tiles_w = src_alignment[ 0 ].size(); - constexpr int repeat_factor = pyramid_scale_factor_prev_curr / tilesize_scale_factor_prev_curr; + dst_alignment[tile_row_i][tile_col_i] = nbr_align_i; + curr_align_i_distance = nbr_align_i_distance; + } + } + } + } + } - // printf("build_upsampled_prev_aligement with scale factor %d, repeat factor %d, tile size factor %d\n", \ - // pyramid_scale_factor_prev_curr, repeat_factor, tilesize_scale_factor_prev_curr ); + } + } - int dst_num_tiles_main_h = src_num_tiles_h * repeat_factor; - int dst_num_tiles_main_w = src_num_tiles_w * repeat_factor; - if ( dst_num_tiles_main_h > num_tiles_h || dst_num_tiles_main_w > num_tiles_w ) - { -#ifdef __ANDROID__ - return; -#else - throw std::runtime_error("current level number of tiles smaller than upsampled tiles\n"); -#endif - } - - // Allocate data for dst_alignment - // NOTE: number of tiles h, number of tiles w might be different from dst_num_tiles_main_h, dst_num_tiles_main_w - // For tiles between num_tile_h and dst_num_tiles_main_h, use (0,0) - dst_alignment.resize( num_tiles_h, std::vector>( num_tiles_w, std::pair(0, 0) ) ); - - // Upsample alignment - #pragma omp parallel for collapse(2) - for ( int row_i = 0; row_i < src_num_tiles_h; row_i++ ) - { - for ( int col_i = 0; col_i < src_num_tiles_w; col_i++ ) - { - // Scale alignment - std::pair align_i = src_alignment[ row_i ][ col_i ]; - align_i.first *= pyramid_scale_factor_prev_curr; - align_i.second *= pyramid_scale_factor_prev_curr; - - // repeat - UNROLL_LOOP( repeat_factor ) - for ( int repeat_row_i = 0; repeat_row_i < repeat_factor; ++repeat_row_i ) - { - int repeat_row_i_offset = row_i * repeat_factor + repeat_row_i; - UNROLL_LOOP( repeat_factor ) - for ( int repeat_col_i = 0; repeat_col_i < repeat_factor; ++repeat_col_i ) - { - int repeat_col_i_offset = col_i * repeat_factor + repeat_col_i; - dst_alignment[ repeat_row_i_offset ][ repeat_col_i_offset ] = align_i; - } - } - } - } - - if ( consider_nbr ) - { - // Copy consurtctor - std::vector>> upsampled_alignment{ dst_alignment }; - - // Distance function - unsigned long long (*distance_func_ptr)(const cv::Mat&, const cv::Mat&, int, int, int, int) = \ - &l1_distance; - - #pragma omp parallel for collapse(2) - for ( int tile_row_i = 0; tile_row_i < num_tiles_h; tile_row_i++ ) - { - for ( int tile_col_i = 0; tile_col_i < num_tiles_w; tile_col_i++ ) - { - const auto& curr_align_i = upsampled_alignment[ tile_row_i ][ tile_col_i ]; - - // Container for nbr alignment pair - std::vector> nbrs_align_i; - - // Consider 4 neighbour's alignment - // Only compute distance if alignment is different - if ( tile_col_i > 0 ) - { - const auto& nbr1_align_i = upsampled_alignment[ tile_row_i + 0 ][ tile_col_i - 1 ]; - if ( curr_align_i != nbr1_align_i ) nbrs_align_i.emplace_back( nbr1_align_i ); - } - - if ( tile_col_i < num_tiles_w - 1 ) - { - const auto& nbr2_align_i = upsampled_alignment[ tile_row_i + 0 ][ tile_col_i + 1 ]; - if ( curr_align_i != nbr2_align_i ) nbrs_align_i.emplace_back( nbr2_align_i ); - } - - if ( tile_row_i > 0 ) - { - const auto& nbr3_align_i = upsampled_alignment[ tile_row_i - 1 ][ tile_col_i + 0 ]; - if ( curr_align_i != nbr3_align_i ) nbrs_align_i.emplace_back( nbr3_align_i ); - } - - if ( tile_row_i < num_tiles_h - 1 ) - { - const auto& nbr4_align_i = upsampled_alignment[ tile_row_i + 1 ][ tile_col_i + 0 ]; - if ( curr_align_i != nbr4_align_i ) nbrs_align_i.emplace_back( nbr4_align_i ); - } - - // If there is a nbr alignment that need to be considered. Compute distance - if ( ! nbrs_align_i.empty() ) - { - int ref_tile_row_start_idx_i = tile_row_i * tile_size / 2; - int ref_tile_col_start_idx_i = tile_col_i * tile_size / 2; - - // curr_align_i's distance - auto curr_align_i_distance = distance_func_ptr( - ref_img, alt_img, \ - ref_tile_row_start_idx_i, \ - ref_tile_col_start_idx_i, \ - ref_tile_row_start_idx_i + curr_align_i.first, \ - ref_tile_col_start_idx_i + curr_align_i.second ); - - for ( const auto& nbr_align_i : nbrs_align_i ) - { - auto nbr_align_i_distance = distance_func_ptr( - ref_img, alt_img, \ - ref_tile_row_start_idx_i, \ - ref_tile_col_start_idx_i, \ - ref_tile_row_start_idx_i + nbr_align_i.first, \ - ref_tile_col_start_idx_i + nbr_align_i.second ); - - if ( nbr_align_i_distance < curr_align_i_distance ) - { - #ifdef NDEBUG - printf("tile [%d, %d] update align, prev align (%d, %d) curr align (%d, %d), prev distance %d curr distance %d\n", \ - tile_row_i, tile_col_i, \ - curr_align_i.first, curr_align_i.second, \ - nbr_align_i.first, nbr_align_i.second, \ - int(curr_align_i_distance), int(nbr_align_i_distance) ); - #endif - - dst_alignment[ tile_row_i ][ tile_col_i ] = nbr_align_i; - curr_align_i_distance = nbr_align_i_distance; - } - } - } - } - } - - } -} - - -// Set tilesize as template argument for better compiler optimization result. -template< typename data_type, typename return_type, int tile_size > -static unsigned long long l1_distance( const cv::Mat& img1, const cv::Mat& img2, \ - int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ - int img2_tile_row_start_idx, int img2_tile_col_start_idx ) -{ - #define CUSTOME_ABS( x ) ( x ) > 0 ? ( x ) : - ( x ) + // Set tilesize as template argument for better compiler optimization result. + template< typename data_type, typename return_type, int tile_size > + static unsigned long long l1_distance(const cv::Mat& img1, const cv::Mat& img2, \ + int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ + int img2_tile_row_start_idx, int img2_tile_col_start_idx) + { +#define CUSTOME_ABS( x ) ( x ) > 0 ? ( x ) : - ( x ) - const data_type* img1_ptr = (const data_type*)img1.data; - const data_type* img2_ptr = (const data_type*)img2.data; + const data_type* img1_ptr = (const data_type*)img1.data; + const data_type* img2_ptr = (const data_type*)img2.data; - int img1_step = img1.step1(); - int img2_step = img2.step1(); + int img1_step = img1.step1(); + int img2_step = img2.step1(); - int img1_width = img1.size().width; - int img1_height = img1.size().height; + int img1_width = img1.size().width; + int img1_height = img1.size().height; - int img2_width = img2.size().width; - int img2_height = img2.size().height; + int img2_width = img2.size().width; + int img2_height = img2.size().height; - // Range check for safety - if ( img1_tile_row_start_idx < 0 || img1_tile_row_start_idx > img1_height - tile_size ) - { + // Range check for safety + if (img1_tile_row_start_idx < 0 || img1_tile_row_start_idx > img1_height - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l1 distance img1_tile_row_start_idx" + std::to_string( img1_tile_row_start_idx ) + \ - " out of valid range (0, " + std::to_string( img1_height - tile_size ) + ")\n" ); + throw std::runtime_error("l1 distance img1_tile_row_start_idx" + std::to_string(img1_tile_row_start_idx) + \ + " out of valid range (0, " + std::to_string(img1_height - tile_size) + ")\n"); #endif - } + } - if ( img1_tile_col_start_idx < 0 || img1_tile_col_start_idx > img1_width - tile_size ) - { + if (img1_tile_col_start_idx < 0 || img1_tile_col_start_idx > img1_width - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l1 distance img1_tile_col_start_idx" + std::to_string( img1_tile_col_start_idx ) + \ - " out of valid range (0, " + std::to_string( img1_width - tile_size ) + ")\n" ); + throw std::runtime_error("l1 distance img1_tile_col_start_idx" + std::to_string(img1_tile_col_start_idx) + \ + " out of valid range (0, " + std::to_string(img1_width - tile_size) + ")\n"); #endif - } + } - if ( img2_tile_row_start_idx < 0 || img2_tile_row_start_idx > img2_height - tile_size ) - { + if (img2_tile_row_start_idx < 0 || img2_tile_row_start_idx > img2_height - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l1 distance img2_tile_row_start_idx out of valid range\n"); + throw std::runtime_error("l1 distance img2_tile_row_start_idx out of valid range\n"); #endif - } + } - if ( img2_tile_col_start_idx < 0 || img2_tile_col_start_idx > img2_width - tile_size ) - { + if (img2_tile_col_start_idx < 0 || img2_tile_col_start_idx > img2_width - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l1 distance img2_tile_col_start_idx out of valid range\n"); + throw std::runtime_error("l1 distance img2_tile_col_start_idx out of valid range\n"); #endif - } + } - return_type sum(0); + return_type sum(0); - UNROLL_LOOP( tile_size ) - for ( int row_i = 0; row_i < tile_size; ++row_i ) - { - const data_type* img1_ptr_row_i = img1_ptr + (img1_tile_row_start_idx + row_i) * img1_step + img1_tile_col_start_idx; - const data_type* img2_ptr_row_i = img2_ptr + (img2_tile_row_start_idx + row_i) * img2_step + img2_tile_col_start_idx; + UNROLL_LOOP(tile_size) + for (int row_i = 0; row_i < tile_size; ++row_i) + { + const data_type* img1_ptr_row_i = img1_ptr + (img1_tile_row_start_idx + row_i) * img1_step + img1_tile_col_start_idx; + const data_type* img2_ptr_row_i = img2_ptr + (img2_tile_row_start_idx + row_i) * img2_step + img2_tile_col_start_idx; - UNROLL_LOOP( tile_size ) - for ( int col_i = 0; col_i < tile_size; ++col_i ) - { - data_type l1 = CUSTOME_ABS( img1_ptr_row_i[ col_i ] - img2_ptr_row_i[ col_i ] ); - sum += l1; - } - } + UNROLL_LOOP(tile_size) + for (int col_i = 0; col_i < tile_size; ++col_i) + { + data_type l1 = CUSTOME_ABS(img1_ptr_row_i[col_i] - img2_ptr_row_i[col_i]); + sum += l1; + } + } - #undef CUSTOME_ABS +#undef CUSTOME_ABS - return sum; -} + return sum; + } -template< typename data_type, typename return_type, int tile_size > -static return_type l2_distance( const cv::Mat& img1, const cv::Mat& img2, \ - int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ - int img2_tile_row_start_idx, int img2_tile_col_start_idx ) -{ - #define CUSTOME_ABS( x ) ( x ) > 0 ? ( x ) : - ( x ) + template< typename data_type, typename return_type, int tile_size > + static return_type l2_distance(const cv::Mat& img1, const cv::Mat& img2, \ + int img1_tile_row_start_idx, int img1_tile_col_start_idx, \ + int img2_tile_row_start_idx, int img2_tile_col_start_idx) + { +#define CUSTOME_ABS( x ) ( x ) > 0 ? ( x ) : - ( x ) - const data_type* img1_ptr = (const data_type*)img1.data; - const data_type* img2_ptr = (const data_type*)img2.data; + const data_type* img1_ptr = (const data_type*)img1.data; + const data_type* img2_ptr = (const data_type*)img2.data; - int img1_step = img1.step1(); - int img2_step = img2.step1(); + int img1_step = img1.step1(); + int img2_step = img2.step1(); - int img1_width = img1.size().width; - int img1_height = img1.size().height; + int img1_width = img1.size().width; + int img1_height = img1.size().height; - int img2_width = img2.size().width; - int img2_height = img2.size().height; + int img2_width = img2.size().width; + int img2_height = img2.size().height; - // Range check for safety - if ( img1_tile_row_start_idx < 0 || img1_tile_row_start_idx > img1_height - tile_size ) - { + // Range check for safety + if (img1_tile_row_start_idx < 0 || img1_tile_row_start_idx > img1_height - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l2 distance img1_tile_row_start_idx" + std::to_string( img1_tile_row_start_idx ) + \ - " out of valid range (0, " + std::to_string( img1_height - tile_size ) + ")\n" ); + throw std::runtime_error("l2 distance img1_tile_row_start_idx" + std::to_string(img1_tile_row_start_idx) + \ + " out of valid range (0, " + std::to_string(img1_height - tile_size) + ")\n"); #endif - } + } - if ( img1_tile_col_start_idx < 0 || img1_tile_col_start_idx > img1_width - tile_size ) - { + if (img1_tile_col_start_idx < 0 || img1_tile_col_start_idx > img1_width - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l2 distance img1_tile_col_start_idx" + std::to_string( img1_tile_col_start_idx ) + \ - " out of valid range (0, " + std::to_string( img1_width - tile_size ) + ")\n" ); + throw std::runtime_error("l2 distance img1_tile_col_start_idx" + std::to_string(img1_tile_col_start_idx) + \ + " out of valid range (0, " + std::to_string(img1_width - tile_size) + ")\n"); #endif - } + } - if ( img2_tile_row_start_idx < 0 || img2_tile_row_start_idx > img2_height - tile_size ) - { + if (img2_tile_row_start_idx < 0 || img2_tile_row_start_idx > img2_height - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l2 distance img2_tile_row_start_idx out of valid range\n"); + throw std::runtime_error("l2 distance img2_tile_row_start_idx out of valid range\n"); #endif - } + } - if ( img2_tile_col_start_idx < 0 || img2_tile_col_start_idx > img2_width - tile_size ) - { + if (img2_tile_col_start_idx < 0 || img2_tile_col_start_idx > img2_width - tile_size) + { #ifdef __ANDROID__ - return 0; + return 0; #else - throw std::runtime_error("l2 distance img2_tile_col_start_idx out of valid range\n"); + throw std::runtime_error("l2 distance img2_tile_col_start_idx out of valid range\n"); #endif - } + } - // printf("Search two tile with ref : \n"); - // print_tile( img1, tile_size, img1_tile_row_start_idx, img1_tile_col_start_idx ); - // printf("Search two tile with alt :\n"); - // print_tile( img2, tile_size, img2_tile_row_start_idx, img2_tile_col_start_idx ); + // printf("Search two tile with ref : \n"); + // print_tile( img1, tile_size, img1_tile_row_start_idx, img1_tile_col_start_idx ); + // printf("Search two tile with alt :\n"); + // print_tile( img2, tile_size, img2_tile_row_start_idx, img2_tile_col_start_idx ); - return_type sum(0); + return_type sum(0); - UNROLL_LOOP( tile_size ) - for ( int row_i = 0; row_i < tile_size; ++row_i ) - { - const data_type* img1_ptr_row_i = img1_ptr + (img1_tile_row_start_idx + row_i) * img1_step + img1_tile_col_start_idx; - const data_type* img2_ptr_row_i = img2_ptr + (img2_tile_row_start_idx + row_i) * img2_step + img2_tile_col_start_idx; + UNROLL_LOOP(tile_size) + for (int row_i = 0; row_i < tile_size; ++row_i) + { + const data_type* img1_ptr_row_i = img1_ptr + (img1_tile_row_start_idx + row_i) * img1_step + img1_tile_col_start_idx; + const data_type* img2_ptr_row_i = img2_ptr + (img2_tile_row_start_idx + row_i) * img2_step + img2_tile_col_start_idx; - UNROLL_LOOP( tile_size ) - for ( int col_i = 0; col_i < tile_size; ++col_i ) - { - data_type l1 = CUSTOME_ABS( img1_ptr_row_i[ col_i ] - img2_ptr_row_i[ col_i ] ); - sum += ( l1 * l1 ); - } - } + UNROLL_LOOP(tile_size) + for (int col_i = 0; col_i < tile_size; ++col_i) + { + data_type l1 = CUSTOME_ABS(img1_ptr_row_i[col_i] - img2_ptr_row_i[col_i]); + sum += (l1 * l1); + } + } - #undef CUSTOME_ABS +#undef CUSTOME_ABS - return sum; -} + return sum; + } -template -static cv::Mat extract_img_tile( const cv::Mat& img, int img_tile_row_start_idx, int img_tile_col_start_idx ) -{ - const T* img_ptr = (const T*)img.data; - int img_width = img.size().width; - int img_height = img.size().height; - int img_step = img.step1(); + template + static cv::Mat extract_img_tile(const cv::Mat& img, int img_tile_row_start_idx, int img_tile_col_start_idx) + { + const T* img_ptr = (const T*)img.data; + int img_width = img.size().width; + int img_height = img.size().height; + int img_step = img.step1(); - if ( img_tile_row_start_idx < 0 || img_tile_row_start_idx > img_height - tile_size ) - { + if (img_tile_row_start_idx < 0 || img_tile_row_start_idx > img_height - tile_size) + { #ifdef __ANDROID__ - return cv::Mat(); + return cv::Mat(); #else - throw std::runtime_error("extract_img_tile img_tile_row_start_idx " + std::to_string( img_tile_row_start_idx ) + \ - " out of valid range (0, " + std::to_string( img_height - tile_size ) + ")\n" ); + throw std::runtime_error("extract_img_tile img_tile_row_start_idx " + std::to_string(img_tile_row_start_idx) + \ + " out of valid range (0, " + std::to_string(img_height - tile_size) + ")\n"); #endif - } + } - if ( img_tile_col_start_idx < 0 || img_tile_col_start_idx > img_width - tile_size ) - { + if (img_tile_col_start_idx < 0 || img_tile_col_start_idx > img_width - tile_size) + { #ifdef __ANDROID__ - return cv::Mat(); + return cv::Mat(); #else - throw std::runtime_error("extract_img_tile img_tile_col_start_idx " + std::to_string( img_tile_col_start_idx ) + \ - " out of valid range (0, " + std::to_string( img_width - tile_size ) + ")\n" ); + throw std::runtime_error("extract_img_tile img_tile_col_start_idx " + std::to_string(img_tile_col_start_idx) + \ + " out of valid range (0, " + std::to_string(img_width - tile_size) + ")\n"); #endif - } - - cv::Mat img_tile( tile_size, tile_size, img.type() ); - T* img_tile_ptr = (T*)img_tile.data; - int img_tile_step = img_tile.step1(); - - UNROLL_LOOP( tile_size ) - for ( int row_i = 0; row_i < tile_size; ++row_i ) - { - const T* img_ptr_row_i = img_ptr + img_step * ( img_tile_row_start_idx + row_i ); - T* img_tile_ptr_row_i = img_tile_ptr + img_tile_step * row_i; - - UNROLL_LOOP( tile_size ) - for ( int col_i = 0; col_i < tile_size; ++col_i ) - { - img_tile_ptr_row_i[ col_i ] = img_ptr_row_i[ img_tile_col_start_idx + col_i ]; - } - } - - return img_tile; -} - - -void align_image_level( \ - const cv::Mat& ref_img, \ - const cv::Mat& alt_img, \ - std::vector>>& prev_aligement, \ - std::vector>>& curr_alignment, \ - int scale_factor_prev_curr, \ - int curr_tile_size, \ - int prev_tile_size, \ - int search_radiou, \ - int distance_type ) -{ - // Every align image level share the same distance function. - // Use function ptr to reduce if else overhead inside for loop - unsigned long long (*distance_func_ptr)(const cv::Mat&, const cv::Mat&, int, int, int, int) = nullptr; - - if ( distance_type == 1 ) // l1 distance - { - if ( curr_tile_size == 8 ) - { - distance_func_ptr = &l1_distance; - } - else if ( curr_tile_size == 16 ) - { - distance_func_ptr = &l1_distance; - } - } - else if ( distance_type == 2 ) // l2 distance - { - if ( curr_tile_size == 8 ) - { - distance_func_ptr = &l2_distance; - } - else if ( curr_tile_size == 16 ) - { - distance_func_ptr = &l2_distance; - } - } - - // Every level share the same upsample function - void (*upsample_alignment_func_ptr)(const std::vector>>&, \ - std::vector>>&, \ - int, int, const cv::Mat&, const cv::Mat&, bool) = nullptr; - if ( scale_factor_prev_curr == 2 ) - { - if ( curr_tile_size / prev_tile_size == 2 ) - { - if ( curr_tile_size == 8 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 2, 8>; - } - else if ( curr_tile_size == 16 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 2, 16>; - } - else - { + } + + cv::Mat img_tile(tile_size, tile_size, img.type()); + T* img_tile_ptr = (T*)img_tile.data; + int img_tile_step = img_tile.step1(); + + UNROLL_LOOP(tile_size) + for (int row_i = 0; row_i < tile_size; ++row_i) + { + const T* img_ptr_row_i = img_ptr + img_step * (img_tile_row_start_idx + row_i); + T* img_tile_ptr_row_i = img_tile_ptr + img_tile_step * row_i; + + UNROLL_LOOP(tile_size) + for (int col_i = 0; col_i < tile_size; ++col_i) + { + img_tile_ptr_row_i[col_i] = img_ptr_row_i[img_tile_col_start_idx + col_i]; + } + } + + return img_tile; + } + + + void align_image_level(\ + const cv::Mat& ref_img, \ + const cv::Mat& alt_img, \ + std::vector>>& prev_aligement, \ + std::vector>>& curr_alignment, \ + int scale_factor_prev_curr, \ + int curr_tile_size, \ + int prev_tile_size, \ + int search_radiou, \ + int distance_type) + { + // Every align image level share the same distance function. + // Use function ptr to reduce if else overhead inside for loop + unsigned long long(*distance_func_ptr)(const cv::Mat&, const cv::Mat&, int, int, int, int) = nullptr; + + if (distance_type == 1) // l1 distance + { + if (curr_tile_size == 8) + { + distance_func_ptr = &l1_distance; + } + else if (curr_tile_size == 16) + { + distance_func_ptr = &l1_distance; + } + } + else if (distance_type == 2) // l2 distance + { + if (curr_tile_size == 8) + { + distance_func_ptr = &l2_distance; + } + else if (curr_tile_size == 16) + { + distance_func_ptr = &l2_distance; + } + } + + // Every level share the same upsample function + void(*upsample_alignment_func_ptr)(const std::vector>>&, \ + std::vector>>&, \ + int, int, const cv::Mat&, const cv::Mat&, bool) = nullptr; + if (scale_factor_prev_curr == 2) + { + if (curr_tile_size / prev_tile_size == 2) + { + if (curr_tile_size == 8) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 2, 8>; + } + else if (curr_tile_size == 16) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 2, 16>; + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); #endif - } - - } - else if ( curr_tile_size / prev_tile_size == 1 ) - { - if ( curr_tile_size == 8 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 1, 8>; - } - else if ( curr_tile_size == 16 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 1, 16>; - } - else - { + } + + } + else if (curr_tile_size / prev_tile_size == 1) + { + if (curr_tile_size == 8) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 1, 8>; + } + else if (curr_tile_size == 16) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<2, 1, 16>; + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); #endif - } - } - else - { + } + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); #endif - } - } - else if ( scale_factor_prev_curr == 4 ) - { - if ( curr_tile_size / prev_tile_size == 2 ) - { - if ( curr_tile_size == 8 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 2, 8>; - } - else if ( curr_tile_size == 16 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 2, 16>; - } - else - { + } + } + else if (scale_factor_prev_curr == 4) + { + if (curr_tile_size / prev_tile_size == 2) + { + if (curr_tile_size == 8) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 2, 8>; + } + else if (curr_tile_size == 16) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 2, 16>; + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); #endif - } - - } - else if ( curr_tile_size / prev_tile_size == 1 ) - { - if ( curr_tile_size == 8 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 1, 8>; - } - else if ( curr_tile_size == 16 ) - { - upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 1, 16>; - } - else - { + } + + } + else if (curr_tile_size / prev_tile_size == 1) + { + if (curr_tile_size == 8) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 1, 8>; + } + else if (curr_tile_size == 16) + { + upsample_alignment_func_ptr = &build_upsampled_prev_aligement<4, 1, 16>; + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); #endif - } - } - else - { + } + } + else + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Something wrong with upsampling function setting\n"); + throw std::runtime_error("Something wrong with upsampling function setting\n"); +#endif + } + } + + + // Function to extract reference image tile for memory cache + cv::Mat(*extract_ref_img_tile)(const cv::Mat&, int, int) = nullptr; + if (curr_tile_size == 8) + { + extract_ref_img_tile = &extract_img_tile; + } + else if (curr_tile_size == 16) + { + extract_ref_img_tile = &extract_img_tile; + } + + // Function to extract search image tile for memory cache + cv::Mat(*extract_alt_img_search)(const cv::Mat&, int, int) = nullptr; + if (curr_tile_size == 8) + { + if (search_radiou == 1) + { + extract_alt_img_search = &extract_img_tile; + } + else if (search_radiou == 4) + { + extract_alt_img_search = &extract_img_tile; + } + } + else if (curr_tile_size == 16) + { + if (search_radiou == 1) + { + extract_alt_img_search = &extract_img_tile; + } + else if (search_radiou == 4) + { + extract_alt_img_search = &extract_img_tile; + } + } + + int num_tiles_h = ref_img.size().height / (curr_tile_size / 2) - 1; + int num_tiles_w = ref_img.size().width / (curr_tile_size / 2) - 1; + + /* Upsample pervious layer alignment */ + std::vector>> upsampled_prev_aligement; + + // Coarsest level + // prev_alignment is invalid / empty, construct alignment as (0,0) + if (prev_tile_size == -1) + { + upsampled_prev_aligement.resize(num_tiles_h, \ + std::vector>(num_tiles_w, std::pair(0, 0))); + } + // Upsample previous level alignment + else + { + upsample_alignment_func_ptr(prev_aligement, upsampled_prev_aligement, \ + num_tiles_h, num_tiles_w, ref_img, alt_img, false); + + // printf("\n!!!!!Upsampled previous alignment\n"); + // for ( int tile_row = 0; tile_row < int(upsampled_prev_aligement.size()); tile_row++ ) + // { + // for ( int tile_col = 0; tile_col < int(upsampled_prev_aligement.at(0).size()); tile_col++ ) + // { + // const auto tile_start = upsampled_prev_aligement.at( tile_row ).at( tile_col ); + // printf("up tile (%d, %d) -> start idx (%d, %d)\n", \ + // tile_row, tile_col, tile_start.first, tile_start.second); + // } + // } + + } + +#ifndef NDEBUG + printf("%s::%s start: \n", __FILE__, __func__); + printf(" scale_factor_prev_curr %d, tile_size %d, prev_tile_size %d, search_radiou %d, distance L%d, \n", \ + scale_factor_prev_curr, curr_tile_size, prev_tile_size, search_radiou, distance_type); + printf(" ref img size h=%d w=%d, alt img size h=%d w=%d, \n", \ + ref_img.size().height, ref_img.size().width, alt_img.size().height, alt_img.size().width); + printf(" num tile h (upsampled) %d, num tile w (upsampled) %d\n", num_tiles_h, num_tiles_w); #endif - } - } - - - // Function to extract reference image tile for memory cache - cv::Mat (*extract_ref_img_tile)(const cv::Mat&, int, int) = nullptr; - if ( curr_tile_size == 8 ) - { - extract_ref_img_tile = &extract_img_tile; - } - else if ( curr_tile_size == 16 ) - { - extract_ref_img_tile = &extract_img_tile; - } - - // Function to extract search image tile for memory cache - cv::Mat (*extract_alt_img_search)(const cv::Mat&, int, int) = nullptr; - if ( curr_tile_size == 8 ) - { - if ( search_radiou == 1 ) - { - extract_alt_img_search = &extract_img_tile; - } - else if ( search_radiou == 4 ) - { - extract_alt_img_search = &extract_img_tile; - } - } - else if ( curr_tile_size == 16 ) - { - if ( search_radiou == 1 ) - { - extract_alt_img_search = &extract_img_tile; - } - else if ( search_radiou == 4 ) - { - extract_alt_img_search = &extract_img_tile; - } - } - - int num_tiles_h = ref_img.size().height / (curr_tile_size / 2) - 1; - int num_tiles_w = ref_img.size().width / (curr_tile_size / 2 ) - 1; - - /* Upsample pervious layer alignment */ - std::vector>> upsampled_prev_aligement; - - // Coarsest level - // prev_alignment is invalid / empty, construct alignment as (0,0) - if ( prev_tile_size == -1 ) - { - upsampled_prev_aligement.resize( num_tiles_h, \ - std::vector>( num_tiles_w, std::pair(0, 0) ) ); - } - // Upsample previous level alignment - else - { - upsample_alignment_func_ptr( prev_aligement, upsampled_prev_aligement, \ - num_tiles_h, num_tiles_w, ref_img, alt_img, false ); - - // printf("\n!!!!!Upsampled previous alignment\n"); - // for ( int tile_row = 0; tile_row < int(upsampled_prev_aligement.size()); tile_row++ ) - // { - // for ( int tile_col = 0; tile_col < int(upsampled_prev_aligement.at(0).size()); tile_col++ ) - // { - // const auto tile_start = upsampled_prev_aligement.at( tile_row ).at( tile_col ); - // printf("up tile (%d, %d) -> start idx (%d, %d)\n", \ - // tile_row, tile_col, tile_start.first, tile_start.second); - // } - // } - - } - - #ifndef NDEBUG - printf("%s::%s start: \n", __FILE__, __func__ ); - printf(" scale_factor_prev_curr %d, tile_size %d, prev_tile_size %d, search_radiou %d, distance L%d, \n", \ - scale_factor_prev_curr, curr_tile_size, prev_tile_size, search_radiou, distance_type ); - printf(" ref img size h=%d w=%d, alt img size h=%d w=%d, \n", \ - ref_img.size().height, ref_img.size().width, alt_img.size().height, alt_img.size().width ); - printf(" num tile h (upsampled) %d, num tile w (upsampled) %d\n", num_tiles_h, num_tiles_w); - #endif - - // allocate memory for current alignmenr - curr_alignment.resize( num_tiles_h, std::vector>( num_tiles_w, std::pair(0, 0) ) ); - - /* Pad alternative image */ - cv::Mat alt_img_pad; - cv::copyMakeBorder( alt_img, \ - alt_img_pad, \ - search_radiou, search_radiou, search_radiou, search_radiou, \ - cv::BORDER_CONSTANT, cv::Scalar( UINT_LEAST16_MAX ) ); - - // printf("Reference image h=%d, w=%d: \n", ref_img.size().height, ref_img.size().width ); - // print_img( ref_img ); - - // printf("Alter image pad h=%d, w=%d: \n", alt_img_pad.size().height, alt_img_pad.size().width ); - // print_img( alt_img_pad ); - - // printf("!! enlarged tile size %d\n", curr_tile_size + 2 * search_radiou ); - - int alt_tile_row_idx_max = alt_img_pad.size().height - ( curr_tile_size + 2 * search_radiou ); - int alt_tile_col_idx_max = alt_img_pad.size().width - ( curr_tile_size + 2 * search_radiou ); - - // Dlete below distance vector, this is for debug only - std::vector> distances( num_tiles_h, std::vector( num_tiles_w, 0 )); - - /* Iterate through all reference tile & compute distance */ - #pragma omp parallel for collapse(2) - for ( int ref_tile_row_i = 0; ref_tile_row_i < num_tiles_h; ref_tile_row_i++ ) - { - for ( int ref_tile_col_i = 0; ref_tile_col_i < num_tiles_w; ref_tile_col_i++ ) - { - // Upper left index of reference tile - int ref_tile_row_start_idx_i = ref_tile_row_i * curr_tile_size / 2; - int ref_tile_col_start_idx_i = ref_tile_col_i * curr_tile_size / 2; - - // printf("\nRef img tile [%d, %d] -> start idx [%d, %d] (row, col)\n", \ - // ref_tile_row_i, ref_tile_col_i, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i ); - // printf("\nRef img tile [%d, %d]\n", ref_tile_row_i, ref_tile_col_i ); - // print_tile( ref_img, curr_tile_size, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i ); - - // Upsampled alignment at this tile - // Alignment are relative displacement in pixel value - int prev_alignment_row_i = upsampled_prev_aligement.at( ref_tile_row_i ).at( ref_tile_col_i ).first; - int prev_alignment_col_i = upsampled_prev_aligement.at( ref_tile_row_i ).at( ref_tile_col_i ).second; - - // Alternative image tile start idx - int alt_tile_row_start_idx_i = ref_tile_row_start_idx_i + prev_alignment_row_i; - int alt_tile_col_start_idx_i = ref_tile_col_start_idx_i + prev_alignment_col_i; - - // Ensure alternative image tile within range - if ( alt_tile_row_start_idx_i < 0 ) - alt_tile_row_start_idx_i = 0; - if ( alt_tile_col_start_idx_i < 0 ) - alt_tile_col_start_idx_i = 0; - if ( alt_tile_row_start_idx_i > alt_tile_row_idx_max ) - { - // int before = alt_tile_row_start_idx_i; - alt_tile_row_start_idx_i = alt_tile_row_idx_max; - // printf("@@ change start x from %d to %d\n", before, alt_tile_row_idx_max); - } - if ( alt_tile_col_start_idx_i > alt_tile_col_idx_max ) - { - // int before = alt_tile_col_start_idx_i; - alt_tile_col_start_idx_i = alt_tile_col_idx_max; - // printf("@@ change start y from %d to %d\n", before, alt_tile_col_idx_max ); - } - - // Explicitly caching reference image tile - cv::Mat ref_img_tile_i = extract_ref_img_tile( ref_img, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i ); - cv::Mat alt_img_search_i = extract_alt_img_search( alt_img_pad, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i ); - - // Because alternative image is padded with search radious. - // Using same coordinate with reference image will automatically considered search radious * 2 - // printf("Alt image tile [%d, %d]-> start idx [%d, %d]\n", \ - // ref_tile_row_i, ref_tile_col_i, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i ); - // printf("\nAlt image tile [%d, %d]\n", ref_tile_row_i, ref_tile_col_i ); - // print_tile( alt_img_pad, curr_tile_size + 2 * search_radiou, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i ); - - // Search based on L1/L2 distance + + // allocate memory for current alignmenr + curr_alignment.resize(num_tiles_h, std::vector>(num_tiles_w, std::pair(0, 0))); + + /* Pad alternative image */ + cv::Mat alt_img_pad; + cv::copyMakeBorder(alt_img, \ + alt_img_pad, \ + search_radiou, search_radiou, search_radiou, search_radiou, \ + cv::BORDER_CONSTANT, cv::Scalar(UINT_LEAST16_MAX)); + + // printf("Reference image h=%d, w=%d: \n", ref_img.size().height, ref_img.size().width ); + // print_img( ref_img ); + + // printf("Alter image pad h=%d, w=%d: \n", alt_img_pad.size().height, alt_img_pad.size().width ); + // print_img( alt_img_pad ); + + // printf("!! enlarged tile size %d\n", curr_tile_size + 2 * search_radiou ); + + int alt_tile_row_idx_max = alt_img_pad.size().height - (curr_tile_size + 2 * search_radiou); + int alt_tile_col_idx_max = alt_img_pad.size().width - (curr_tile_size + 2 * search_radiou); + + // Dlete below distance vector, this is for debug only + std::vector> distances(num_tiles_h, std::vector(num_tiles_w, 0)); + + /* Iterate through all reference tile & compute distance */ +#pragma omp parallel for collapse(2) + for (int ref_tile_row_i = 0; ref_tile_row_i < num_tiles_h; ref_tile_row_i++) + { + for (int ref_tile_col_i = 0; ref_tile_col_i < num_tiles_w; ref_tile_col_i++) + { + // Upper left index of reference tile + int ref_tile_row_start_idx_i = ref_tile_row_i * curr_tile_size / 2; + int ref_tile_col_start_idx_i = ref_tile_col_i * curr_tile_size / 2; + + // printf("\nRef img tile [%d, %d] -> start idx [%d, %d] (row, col)\n", \ + // ref_tile_row_i, ref_tile_col_i, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i ); + // printf("\nRef img tile [%d, %d]\n", ref_tile_row_i, ref_tile_col_i ); + // print_tile( ref_img, curr_tile_size, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i ); + + // Upsampled alignment at this tile + // Alignment are relative displacement in pixel value + int prev_alignment_row_i = upsampled_prev_aligement.at(ref_tile_row_i).at(ref_tile_col_i).first; + int prev_alignment_col_i = upsampled_prev_aligement.at(ref_tile_row_i).at(ref_tile_col_i).second; + + // Alternative image tile start idx + int alt_tile_row_start_idx_i = ref_tile_row_start_idx_i + prev_alignment_row_i; + int alt_tile_col_start_idx_i = ref_tile_col_start_idx_i + prev_alignment_col_i; + + // Ensure alternative image tile within range + if (alt_tile_row_start_idx_i < 0) + alt_tile_row_start_idx_i = 0; + if (alt_tile_col_start_idx_i < 0) + alt_tile_col_start_idx_i = 0; + if (alt_tile_row_start_idx_i > alt_tile_row_idx_max) + { + // int before = alt_tile_row_start_idx_i; + alt_tile_row_start_idx_i = alt_tile_row_idx_max; + // printf("@@ change start x from %d to %d\n", before, alt_tile_row_idx_max); + } + if (alt_tile_col_start_idx_i > alt_tile_col_idx_max) + { + // int before = alt_tile_col_start_idx_i; + alt_tile_col_start_idx_i = alt_tile_col_idx_max; + // printf("@@ change start y from %d to %d\n", before, alt_tile_col_idx_max ); + } + + // Explicitly caching reference image tile + cv::Mat ref_img_tile_i = extract_ref_img_tile(ref_img, ref_tile_row_start_idx_i, ref_tile_col_start_idx_i); + cv::Mat alt_img_search_i = extract_alt_img_search(alt_img_pad, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i); + + // Because alternative image is padded with search radious. + // Using same coordinate with reference image will automatically considered search radious * 2 + // printf("Alt image tile [%d, %d]-> start idx [%d, %d]\n", \ + // ref_tile_row_i, ref_tile_col_i, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i ); + // printf("\nAlt image tile [%d, %d]\n", ref_tile_row_i, ref_tile_col_i ); + // print_tile( alt_img_pad, curr_tile_size + 2 * search_radiou, alt_tile_row_start_idx_i, alt_tile_col_start_idx_i ); + + // Search based on L1/L2 distance #ifdef _WIN32 - unsigned long long min_distance_i = ULLONG_MAX; + unsigned long long min_distance_i = ULLONG_MAX; #else - unsigned long long min_distance_i = ULONG_LONG_MAX; + unsigned long long min_distance_i = ULONG_LONG_MAX; #endif - int min_distance_row_i = -1; - int min_distance_col_i = -1; - for ( int search_row_j = 0; search_row_j < ( search_radiou * 2 + 1 ); search_row_j++ ) - { - for ( int search_col_j = 0; search_col_j < ( search_radiou * 2 + 1 ); search_col_j++ ) - { - // printf("\n--->tile at [%d, %d] search (%d, %d)\n", \ - // ref_tile_row_i, ref_tile_col_i, search_row_j - search_radiou, search_col_j - search_radiou ); - - // unsigned long long distance_j = distance_func_ptr( ref_img, alt_img_pad, \ - // ref_tile_row_start_idx_i, ref_tile_col_start_idx_i, \ - // alt_tile_row_start_idx_i + search_row_j, alt_tile_col_start_idx_i + search_col_j ); - - // unsigned long long distance_j = distance_func_ptr( ref_img_tile_i, alt_img_pad, \ - // 0, 0, \ - // alt_tile_row_start_idx_i + search_row_j, alt_tile_col_start_idx_i + search_col_j ); - - unsigned long long distance_j = distance_func_ptr( ref_img_tile_i, alt_img_search_i, \ - 0, 0, \ - search_row_j, search_col_j ); - - // printf("<---tile at [%d, %d] search (%d, %d), new dis %llu, old dis %llu\n", \ - // ref_tile_row_i, ref_tile_col_i, search_row_j - search_radiou, search_col_j - search_radiou, distance_j, min_distance_i ); - - // If this is smaller distance - if ( distance_j < min_distance_i ) - { - min_distance_i = distance_j; - min_distance_col_i = search_col_j; - min_distance_row_i = search_row_j; - } - - // If same value, choose the one closer to the original tile location - if ( distance_j == min_distance_i && min_distance_row_i != -1 && min_distance_col_i != -1 ) - { - int prev_distance_row_2_ref = min_distance_row_i - search_radiou; - int prev_distance_col_2_ref = min_distance_col_i - search_radiou; - int curr_distance_row_2_ref = search_row_j - search_radiou; - int curr_distance_col_2_ref = search_col_j - search_radiou; - - int prev_distance_2_ref_sqr = prev_distance_row_2_ref * prev_distance_row_2_ref + prev_distance_col_2_ref * prev_distance_col_2_ref; - int curr_distance_2_ref_sqr = curr_distance_row_2_ref * curr_distance_row_2_ref + curr_distance_col_2_ref * curr_distance_col_2_ref; - - // previous min distance idx is farther away from ref tile start location - if ( prev_distance_2_ref_sqr > curr_distance_2_ref_sqr ) - { - // printf("@@@ Same distance %d, choose closer one (%d, %d) instead of (%d, %d)\n", \ - // distance_j, search_row_j, search_col_j, min_distance_row_i, min_distance_col_i); - min_distance_col_i = search_col_j; - min_distance_row_i = search_row_j; - } - } - } - } - - // printf("tile at (%d, %d) alignment (%d, %d)\n", \ - // ref_tile_row_i, ref_tile_col_i, min_distance_row_i, min_distance_col_i ); - - int alignment_row_i = prev_alignment_row_i + min_distance_row_i - search_radiou; - int alignment_col_i = prev_alignment_col_i + min_distance_col_i - search_radiou; - - std::pair alignment_i( alignment_row_i, alignment_col_i ); - - // Add min_distance_i's corresbonding idx as min - curr_alignment.at( ref_tile_row_i ).at( ref_tile_col_i ) = alignment_i; - distances.at( ref_tile_row_i ).at( ref_tile_col_i ) = min_distance_i; - } - } - - // printf("\n!!!!!Min distance for each tile \n"); - // for ( int tile_row = 0; tile_row < num_tiles_h; tile_row++ ) - // { - // for ( int tile_col = 0; tile_col < num_tiles_w; ++tile_col ) - // { - // printf("tile (%d, %d) distance %u\n", \ - // tile_row, tile_col, distances.at( tile_row).at(tile_col ) ); - // } - // } - - // printf("\n!!!!!Alignment at current level\n"); - // for ( int tile_row = 0; tile_row < num_tiles_h; tile_row++ ) - // { - // for ( int tile_col = 0; tile_col < num_tiles_w; tile_col++ ) - // { - // const auto tile_start = curr_alignment.at( tile_row ).at( tile_col ); - // printf("tile (%d, %d) -> start idx (%d, %d)\n", \ - // tile_row, tile_col, tile_start.first, tile_start.second); - // } - // } - -} - - - -void align::process( const hdrplus::burst& burst_images, \ - std::vector>>>& images_alignment ) -{ - #ifndef NDEBUG - printf("%s::%s align::process start\n", __FILE__, __func__ ); fflush(stdout); - #endif - - images_alignment.clear(); - images_alignment.resize( burst_images.num_images ); - - // image pyramid per image, per pyramid level - std::vector> per_grayimg_pyramid; - - // printf("!!!!! ref bayer padded\n"); - // print_img( burst_images.bayer_images_pad.at( burst_images.reference_image_idx) ); - // exit(1); - - // printf("!!!!! ref gray padded\n"); - // print_img( burst_images.grayscale_images_pad.at( burst_images.reference_image_idx) ); - // exit(1); - - per_grayimg_pyramid.resize( burst_images.num_images ); - - #pragma omp parallel for - for ( int img_idx = 0; img_idx < burst_images.num_images; ++img_idx ) - { - // per_grayimg_pyramid[ img_idx ][ 0 ] is the original image - // per_grayimg_pyramid[ img_idx ][ 3 ] is the coarsest image - build_per_grayimg_pyramid( per_grayimg_pyramid.at( img_idx ), \ - burst_images.grayscale_images_pad.at( img_idx ), \ - this->inv_scale_factors ); - } - - // #ifndef NDEBUG - // printf("%s::%s build image pyramid of size : ", __FILE__, __func__ ); - // for ( int level_i = 0; level_i < num_levels; ++level_i ) - // { - // printf("(%d, %d) ", per_grayimg_pyramid[ 0 ][ level_i ].size().height, - // per_grayimg_pyramid[ 0 ][ level_i ].size().width ); - // } - // printf("\n"); fflush(stdout); - // #endif - - // print image pyramid - // for ( int level_i; level_i < num_levels; ++level_i ) - // { - // printf("\n\n!!!!! ref gray pyramid level %d img : \n" , level_i ); - // print_img( per_grayimg_pyramid[ burst_images.reference_image_idx ][ level_i ] ); - // } - // exit(-1); - - // Align every image - const std::vector& ref_grayimg_pyramid = per_grayimg_pyramid[ burst_images.reference_image_idx ]; - std::vector>> curr_alignment; - std::vector>> prev_alignment; - for ( int img_idx = 0; img_idx < burst_images.num_images; ++img_idx ) - { - // Do not align with reference image - if ( img_idx == burst_images.reference_image_idx ) - continue; - - const std::vector& alt_grayimg_pyramid = per_grayimg_pyramid[ img_idx ]; - - // Align every level from coarse to grain - // level 0 : finest level, the original image - // level 3 : coarsest level - curr_alignment.clear(); - prev_alignment.clear(); - for ( int level_i = num_levels - 1; level_i >= 0; level_i-- ) // 3,2,1,0 - { - // make curr alignment as previous alignment - prev_alignment.swap( curr_alignment ); - curr_alignment.clear(); - - // printf("\n\n########################align level %d\n", level_i ); - align_image_level( - ref_grayimg_pyramid[ level_i ], // reference image at current level - alt_grayimg_pyramid[ level_i ], // alternative image at current level - prev_alignment, // previous layer alignment - curr_alignment, // current layer alignment - ( level_i == ( num_levels - 1 ) ? -1 : inv_scale_factors[ level_i + 1 ] ), // scale factor between previous layer and current layer. -1 if current layer is the coarsest layer, [-1, 4, 4, 2] - grayimg_tile_sizes[ level_i ], // current level tile size - ( level_i == ( num_levels - 1 ) ? -1 : grayimg_tile_sizes[ level_i + 1 ] ), // previous level tile size - grayimg_search_radious[ level_i ], // search radious - distances[ level_i ] ); // L1/L2 distance - - // printf("@@@Alignment at level %d is h=%d, w=%d", level_i, curr_alignment.size(), curr_alignment.at(0).size() ); - - - } // for pyramid level - - // Alignment at grayscale image - images_alignment.at( img_idx ).swap( curr_alignment ); - - // printf("\n!!!!!Alternative Image Alignment\n"); - // for ( int tile_row = 0; tile_row < images_alignment.at( img_idx ).size(); tile_row++ ) - // { - // for ( int tile_col = 0; tile_col < images_alignment.at( img_idx ).at(0).size(); tile_col++ ) - // { - // const auto tile_start = images_alignment.at( img_idx ).at( tile_row ).at( tile_col ); - // printf("tile (%d, %d) -> start idx (%d, %d)\n", \ - // tile_row, tile_col, tile_start.first, tile_start.second); - // } - // } - - } // for alternative image - - per_grayimg_pyramid.clear(); - -} + int min_distance_row_i = -1; + int min_distance_col_i = -1; + for (int search_row_j = 0; search_row_j < (search_radiou * 2 + 1); search_row_j++) + { + for (int search_col_j = 0; search_col_j < (search_radiou * 2 + 1); search_col_j++) + { + // printf("\n--->tile at [%d, %d] search (%d, %d)\n", \ + // ref_tile_row_i, ref_tile_col_i, search_row_j - search_radiou, search_col_j - search_radiou ); + + // unsigned long long distance_j = distance_func_ptr( ref_img, alt_img_pad, \ + // ref_tile_row_start_idx_i, ref_tile_col_start_idx_i, \ + // alt_tile_row_start_idx_i + search_row_j, alt_tile_col_start_idx_i + search_col_j ); + + // unsigned long long distance_j = distance_func_ptr( ref_img_tile_i, alt_img_pad, \ + // 0, 0, \ + // alt_tile_row_start_idx_i + search_row_j, alt_tile_col_start_idx_i + search_col_j ); + + unsigned long long distance_j = distance_func_ptr(ref_img_tile_i, alt_img_search_i, \ + 0, 0, \ + search_row_j, search_col_j); + + // printf("<---tile at [%d, %d] search (%d, %d), new dis %llu, old dis %llu\n", \ + // ref_tile_row_i, ref_tile_col_i, search_row_j - search_radiou, search_col_j - search_radiou, distance_j, min_distance_i ); + + // If this is smaller distance + if (distance_j < min_distance_i) + { + min_distance_i = distance_j; + min_distance_col_i = search_col_j; + min_distance_row_i = search_row_j; + } + + // If same value, choose the one closer to the original tile location + if (distance_j == min_distance_i && min_distance_row_i != -1 && min_distance_col_i != -1) + { + int prev_distance_row_2_ref = min_distance_row_i - search_radiou; + int prev_distance_col_2_ref = min_distance_col_i - search_radiou; + int curr_distance_row_2_ref = search_row_j - search_radiou; + int curr_distance_col_2_ref = search_col_j - search_radiou; + + int prev_distance_2_ref_sqr = prev_distance_row_2_ref * prev_distance_row_2_ref + prev_distance_col_2_ref * prev_distance_col_2_ref; + int curr_distance_2_ref_sqr = curr_distance_row_2_ref * curr_distance_row_2_ref + curr_distance_col_2_ref * curr_distance_col_2_ref; + + // previous min distance idx is farther away from ref tile start location + if (prev_distance_2_ref_sqr > curr_distance_2_ref_sqr) + { + // printf("@@@ Same distance %d, choose closer one (%d, %d) instead of (%d, %d)\n", \ + // distance_j, search_row_j, search_col_j, min_distance_row_i, min_distance_col_i); + min_distance_col_i = search_col_j; + min_distance_row_i = search_row_j; + } + } + } + } + + // printf("tile at (%d, %d) alignment (%d, %d)\n", \ + // ref_tile_row_i, ref_tile_col_i, min_distance_row_i, min_distance_col_i ); + + int alignment_row_i = prev_alignment_row_i + min_distance_row_i - search_radiou; + int alignment_col_i = prev_alignment_col_i + min_distance_col_i - search_radiou; + + std::pair alignment_i(alignment_row_i, alignment_col_i); + + // Add min_distance_i's corresbonding idx as min + curr_alignment.at(ref_tile_row_i).at(ref_tile_col_i) = alignment_i; + distances.at(ref_tile_row_i).at(ref_tile_col_i) = min_distance_i; + } + } + + // printf("\n!!!!!Min distance for each tile \n"); + // for ( int tile_row = 0; tile_row < num_tiles_h; tile_row++ ) + // { + // for ( int tile_col = 0; tile_col < num_tiles_w; ++tile_col ) + // { + // printf("tile (%d, %d) distance %u\n", \ + // tile_row, tile_col, distances.at( tile_row).at(tile_col ) ); + // } + // } + + // printf("\n!!!!!Alignment at current level\n"); + // for ( int tile_row = 0; tile_row < num_tiles_h; tile_row++ ) + // { + // for ( int tile_col = 0; tile_col < num_tiles_w; tile_col++ ) + // { + // const auto tile_start = curr_alignment.at( tile_row ).at( tile_col ); + // printf("tile (%d, %d) -> start idx (%d, %d)\n", \ + // tile_row, tile_col, tile_start.first, tile_start.second); + // } + // } + + } + + + + void align::process(const hdrplus::burst& burst_images, \ + std::vector>>>& images_alignment) + { +#ifndef NDEBUG + printf("%s::%s align::process start\n", __FILE__, __func__); fflush(stdout); +#endif + + images_alignment.clear(); + images_alignment.resize(burst_images.num_images); + + // image pyramid per image, per pyramid level + std::vector> per_grayimg_pyramid; + + // printf("!!!!! ref bayer padded\n"); + // print_img( burst_images.bayer_images_pad.at( burst_images.reference_image_idx) ); + // exit(1); + + // printf("!!!!! ref gray padded\n"); + // print_img( burst_images.grayscale_images_pad.at( burst_images.reference_image_idx) ); + // exit(1); + + per_grayimg_pyramid.resize(burst_images.num_images); + +#pragma omp parallel for + for (int img_idx = 0; img_idx < burst_images.num_images; ++img_idx) + { + // per_grayimg_pyramid[ img_idx ][ 0 ] is the original image + // per_grayimg_pyramid[ img_idx ][ 3 ] is the coarsest image + build_per_grayimg_pyramid(per_grayimg_pyramid.at(img_idx), \ + burst_images.grayscale_images_pad.at(img_idx), \ + this->inv_scale_factors); + } + + // #ifndef NDEBUG + // printf("%s::%s build image pyramid of size : ", __FILE__, __func__ ); + // for ( int level_i = 0; level_i < num_levels; ++level_i ) + // { + // printf("(%d, %d) ", per_grayimg_pyramid[ 0 ][ level_i ].size().height, + // per_grayimg_pyramid[ 0 ][ level_i ].size().width ); + // } + // printf("\n"); fflush(stdout); + // #endif + + // print image pyramid + // for ( int level_i; level_i < num_levels; ++level_i ) + // { + // printf("\n\n!!!!! ref gray pyramid level %d img : \n" , level_i ); + // print_img( per_grayimg_pyramid[ burst_images.reference_image_idx ][ level_i ] ); + // } + // exit(-1); + + // Align every image + const std::vector& ref_grayimg_pyramid = per_grayimg_pyramid[burst_images.reference_image_idx]; + std::vector>> curr_alignment; + std::vector>> prev_alignment; + for (int img_idx = 0; img_idx < burst_images.num_images; ++img_idx) + { + // Do not align with reference image + if (img_idx == burst_images.reference_image_idx) + continue; + + const std::vector& alt_grayimg_pyramid = per_grayimg_pyramid[img_idx]; + + // Align every level from coarse to grain + // level 0 : finest level, the original image + // level 3 : coarsest level + curr_alignment.clear(); + prev_alignment.clear(); + for (int level_i = num_levels - 1; level_i >= 0; level_i--) // 3,2,1,0 + { + // make curr alignment as previous alignment + prev_alignment.swap(curr_alignment); + curr_alignment.clear(); + + // printf("\n\n########################align level %d\n", level_i ); + align_image_level( + ref_grayimg_pyramid[level_i], // reference image at current level + alt_grayimg_pyramid[level_i], // alternative image at current level + prev_alignment, // previous layer alignment + curr_alignment, // current layer alignment + (level_i == (num_levels - 1) ? -1 : inv_scale_factors[level_i + 1]), // scale factor between previous layer and current layer. -1 if current layer is the coarsest layer, [-1, 4, 4, 2] + grayimg_tile_sizes[level_i], // current level tile size + (level_i == (num_levels - 1) ? -1 : grayimg_tile_sizes[level_i + 1]), // previous level tile size + grayimg_search_radious[level_i], // search radious + distances[level_i]); // L1/L2 distance + + // printf("@@@Alignment at level %d is h=%d, w=%d", level_i, curr_alignment.size(), curr_alignment.at(0).size() ); + + + } // for pyramid level + + // Alignment at grayscale image + images_alignment.at(img_idx).swap(curr_alignment); + + // printf("\n!!!!!Alternative Image Alignment\n"); + // for ( int tile_row = 0; tile_row < images_alignment.at( img_idx ).size(); tile_row++ ) + // { + // for ( int tile_col = 0; tile_col < images_alignment.at( img_idx ).at(0).size(); tile_col++ ) + // { + // const auto tile_start = images_alignment.at( img_idx ).at( tile_row ).at( tile_col ); + // printf("tile (%d, %d) -> start idx (%d, %d)\n", \ + // tile_row, tile_col, tile_start.first, tile_start.second); + // } + // } + + } // for alternative image + + per_grayimg_pyramid.clear(); + + } } // namespace hdrplus diff --git a/src/bayer_image.cpp b/src/bayer_image.cpp index a65811a..dd6b182 100644 --- a/src/bayer_image.cpp +++ b/src/bayer_image.cpp @@ -12,227 +12,227 @@ namespace hdrplus { -bayer_image::bayer_image( const std::string& bayer_image_path ) -{ - libraw_processor = std::make_shared(); - - // Open RAW image file - int return_code; - if ( ( return_code = libraw_processor->open_file( bayer_image_path.c_str() ) ) != LIBRAW_SUCCESS ) - { - libraw_processor->recycle(); + bayer_image::bayer_image(const std::string& bayer_image_path) + { + libraw_processor = std::make_shared(); + + // Open RAW image file + int return_code; + if ((return_code = libraw_processor->open_file(bayer_image_path.c_str())) != LIBRAW_SUCCESS) + { + libraw_processor->recycle(); #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror( return_code )); + throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror(return_code)); #endif - } + } - // Unpack the raw image - if ( ( return_code = libraw_processor->unpack() ) != LIBRAW_SUCCESS ) - { + // Unpack the raw image + if ((return_code = libraw_processor->unpack()) != LIBRAW_SUCCESS) + { #ifdef __ANDROID__ - return; + return; #else - throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror( return_code )); + throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror(return_code)); #endif - } - - // Get image basic info - width = int( libraw_processor->imgdata.rawdata.sizes.raw_width ); - height = int( libraw_processor->imgdata.rawdata.sizes.raw_height ); - - // Read exif tags - Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(bayer_image_path); - assert(image.get() != 0); - image->readMetadata(); - Exiv2::ExifData &exifData = image->exifData(); - if (exifData.empty()) { - std::string error(bayer_image_path); - error += ": No Exif data found in the file"; - std::cout << error << std::endl; - } - - white_level = exifData["Exif.Image.WhiteLevel"].toLong(); - black_level_per_channel.resize( 4 ); - black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); - black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); - black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); - black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); - iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); - - // Create CV mat - // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ - // https://www.libraw.org/node/2141 - raw_image = cv::Mat( height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image ).clone(); // changed the order of width and height - - // 2x2 box filter - grayscale_image = box_filter_kxk( raw_image ); - - #ifndef NDEBUG - printf("%s::%s read bayer image %s with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ - __FILE__, __func__, bayer_image_path.c_str(), width, height, iso, white_level, \ - black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3] ); - fflush( stdout ); - #endif -} + } + + // Get image basic info + width = int(libraw_processor->imgdata.rawdata.sizes.raw_width); + height = int(libraw_processor->imgdata.rawdata.sizes.raw_height); + + // Read exif tags + Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(bayer_image_path); + assert(image.get() != 0); + image->readMetadata(); + Exiv2::ExifData &exifData = image->exifData(); + if (exifData.empty()) { + std::string error(bayer_image_path); + error += ": No Exif data found in the file"; + std::cout << error << std::endl; + } + + white_level = exifData["Exif.Image.WhiteLevel"].toLong(); + black_level_per_channel.resize(4); + black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); + black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); + black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); + black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); + iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); + + // Create CV mat + // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ + // https://www.libraw.org/node/2141 + raw_image = cv::Mat(height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image).clone(); // changed the order of width and height + + // 2x2 box filter + grayscale_image = box_filter_kxk(raw_image); -bayer_image::bayer_image( const std::vector& bayer_image_content ) -{ - libraw_processor = std::make_shared(); +#ifndef NDEBUG + printf("%s::%s read bayer image %s with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ + __FILE__, __func__, bayer_image_path.c_str(), width, height, iso, white_level, \ + black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3]); + fflush(stdout); +#endif + } - // Open RAW image file - int return_code; - if ( ( return_code = libraw_processor->open_buffer( (void *)(&bayer_image_content[0]), bayer_image_content.size() ) ) != LIBRAW_SUCCESS ) - { - libraw_processor->recycle(); + bayer_image::bayer_image(const std::vector& bayer_image_content) + { + libraw_processor = std::make_shared(); + + // Open RAW image file + int return_code; + if ((return_code = libraw_processor->open_buffer((void *)(&bayer_image_content[0]), bayer_image_content.size())) != LIBRAW_SUCCESS) + { + libraw_processor->recycle(); #ifdef __ANDROID__ || _WIN32 - return; + return; #else - // throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror( return_code )); - return; + // throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror( return_code )); + return; #endif - } + } - // Unpack the raw image - if ( ( return_code = libraw_processor->unpack() ) != LIBRAW_SUCCESS ) - { + // Unpack the raw image + if ((return_code = libraw_processor->unpack()) != LIBRAW_SUCCESS) + { #ifdef __ANDROID__ || _WIN32 - return; + return; #else - // throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror( return_code )); - return; + // throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror( return_code )); + return; #endif - } - - // Get image basic info - width = int( libraw_processor->imgdata.rawdata.sizes.raw_width ); - height = int( libraw_processor->imgdata.rawdata.sizes.raw_height ); - - // Read exif tags - Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(&bayer_image_content[0], bayer_image_content.size()); - assert(image.get() != 0); - image->readMetadata(); - Exiv2::ExifData &exifData = image->exifData(); - if (exifData.empty()) { - std::string error = "No Exif data found in the file"; - std::cout << error << std::endl; - } - - white_level = exifData["Exif.Image.WhiteLevel"].toLong(); - black_level_per_channel.resize( 4 ); - black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); - black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); - black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); - black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); - iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); - - // Create CV mat - // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ - // https://www.libraw.org/node/2141 - raw_image = cv::Mat( height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image ).clone(); // changed the order of width and height - - // 2x2 box filter - grayscale_image = box_filter_kxk( raw_image ); + } + + // Get image basic info + width = int(libraw_processor->imgdata.rawdata.sizes.raw_width); + height = int(libraw_processor->imgdata.rawdata.sizes.raw_height); + + // Read exif tags + Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(&bayer_image_content[0], bayer_image_content.size()); + assert(image.get() != 0); + image->readMetadata(); + Exiv2::ExifData &exifData = image->exifData(); + if (exifData.empty()) { + std::string error = "No Exif data found in the file"; + std::cout << error << std::endl; + } + + white_level = exifData["Exif.Image.WhiteLevel"].toLong(); + black_level_per_channel.resize(4); + black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); + black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); + black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); + black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); + iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); + + // Create CV mat + // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ + // https://www.libraw.org/node/2141 + raw_image = cv::Mat(height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image).clone(); // changed the order of width and height + + // 2x2 box filter + grayscale_image = box_filter_kxk(raw_image); #ifndef NDEBUG - printf("%s::%s read bayer image with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ - __FILE__, __func__, width, height, iso, white_level, \ - black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3] ); - fflush( stdout ); + printf("%s::%s read bayer image with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ + __FILE__, __func__, width, height, iso, white_level, \ + black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3]); + fflush(stdout); #endif -} - - bayer_image::bayer_image( std::shared_ptr bayer_image_file ) - { - libraw_processor = std::make_shared(); - - // Open RAW image file - int return_code; - { - std::vector& fileData = bayer_image_file->content; - if ( ( return_code = libraw_processor->open_buffer( (void *)(&fileData[0]), fileData.size() ) ) != LIBRAW_SUCCESS ) - { - libraw_processor->recycle(); + } + + bayer_image::bayer_image(std::shared_ptr bayer_image_file) + { + libraw_processor = std::make_shared(); + + // Open RAW image file + int return_code; + { + std::vector& fileData = bayer_image_file->content; + if ((return_code = libraw_processor->open_buffer((void *)(&fileData[0]), fileData.size())) != LIBRAW_SUCCESS) + { + libraw_processor->recycle(); #ifdef __ANDROID__ - return; + return; #else return; - // throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror( return_code )); + // throw std::runtime_error("Error opening file " + bayer_image_path + " " + libraw_strerror( return_code )); #endif - } - } + } + } - // Unpack the raw image - if ( ( return_code = libraw_processor->unpack() ) != LIBRAW_SUCCESS ) - { + // Unpack the raw image + if ((return_code = libraw_processor->unpack()) != LIBRAW_SUCCESS) + { #ifdef __ANDROID__ - return; + return; #else return; - // throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror( return_code )); + // throw std::runtime_error("Error unpack file " + bayer_image_path + " " + libraw_strerror( return_code )); #endif - } - - // Get image basic info - width = int( libraw_processor->imgdata.rawdata.sizes.raw_width ); - height = int( libraw_processor->imgdata.rawdata.sizes.raw_height ); - - // Read exif tags - Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(&bayer_image_file->content[0], bayer_image_file->content.size()); - assert(image.get() != 0); - image->readMetadata(); - Exiv2::ExifData &exifData = image->exifData(); - if (exifData.empty()) { - std::string error = "No Exif data found in the file"; - std::cout << error << std::endl; - } - - white_level = exifData["Exif.Image.WhiteLevel"].toLong(); - black_level_per_channel.resize( 4 ); - black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); - black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); - black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); - black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); - iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); - - // Create CV mat - // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ - // https://www.libraw.org/node/2141 - raw_image = cv::Mat( height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image ).clone(); // changed the order of width and height - - // 2x2 box filter - grayscale_image = box_filter_kxk( raw_image ); + } + + // Get image basic info + width = int(libraw_processor->imgdata.rawdata.sizes.raw_width); + height = int(libraw_processor->imgdata.rawdata.sizes.raw_height); + + // Read exif tags + Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(&bayer_image_file->content[0], bayer_image_file->content.size()); + assert(image.get() != 0); + image->readMetadata(); + Exiv2::ExifData &exifData = image->exifData(); + if (exifData.empty()) { + std::string error = "No Exif data found in the file"; + std::cout << error << std::endl; + } + + white_level = exifData["Exif.Image.WhiteLevel"].toLong(); + black_level_per_channel.resize(4); + black_level_per_channel.at(0) = exifData["Exif.Image.BlackLevel"].toLong(0); + black_level_per_channel.at(1) = exifData["Exif.Image.BlackLevel"].toLong(1); + black_level_per_channel.at(2) = exifData["Exif.Image.BlackLevel"].toLong(2); + black_level_per_channel.at(3) = exifData["Exif.Image.BlackLevel"].toLong(3); + iso = exifData["Exif.Image.ISOSpeedRatings"].toLong(); + + // Create CV mat + // https://answers.opencv.org/question/105972/de-bayering-a-cr2-image/ + // https://www.libraw.org/node/2141 + raw_image = cv::Mat(height, width, CV_16U, libraw_processor->imgdata.rawdata.raw_image).clone(); // changed the order of width and height + + // 2x2 box filter + grayscale_image = box_filter_kxk(raw_image); #ifndef NDEBUG - printf("%s::%s read bayer image with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ - __FILE__, __func__, width, height, iso, white_level, \ - black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3] ); - fflush( stdout ); + printf("%s::%s read bayer image with\n width %zu\n height %zu\n iso %.3f\n white level %d\n black level %d %d %d %d\n", \ + __FILE__, __func__, width, height, iso, white_level, \ + black_level_per_channel[0], black_level_per_channel[1], black_level_per_channel[2], black_level_per_channel[3]); + fflush(stdout); #endif - } + } -std::pair bayer_image::get_noise_params() const -{ - // Set ISO to 100 if not positive - double iso_ = iso <= 0 ? 100 : iso; + std::pair bayer_image::get_noise_params() const + { + // Set ISO to 100 if not positive + double iso_ = iso <= 0 ? 100 : iso; - // Calculate shot noise and read noise parameters w.r.t ISO 100 - double lambda_shot_p = iso_ / 100.0f * baseline_lambda_shot; - double lambda_read_p = (iso_ / 100.0f) * (iso_ / 100.0f) * baseline_lambda_read; + // Calculate shot noise and read noise parameters w.r.t ISO 100 + double lambda_shot_p = iso_ / 100.0f * baseline_lambda_shot; + double lambda_read_p = (iso_ / 100.0f) * (iso_ / 100.0f) * baseline_lambda_read; - double black_level = (black_level_per_channel[0] + \ - black_level_per_channel[1] + \ - black_level_per_channel[2] + \ - black_level_per_channel[3]) / 4.0; + double black_level = (black_level_per_channel[0] + \ + black_level_per_channel[1] + \ + black_level_per_channel[2] + \ + black_level_per_channel[3]) / 4.0; - // Rescale shot and read noise to normal range - double lambda_shot = lambda_shot_p * (white_level - black_level); - double lambda_read = lambda_read_p * (white_level - black_level) * (white_level - black_level); + // Rescale shot and read noise to normal range + double lambda_shot = lambda_shot_p * (white_level - black_level); + double lambda_read = lambda_read_p * (white_level - black_level) * (white_level - black_level); - // return pair - return std::make_pair(lambda_shot, lambda_read); -} + // return pair + return std::make_pair(lambda_shot, lambda_read); + } } diff --git a/src/finish.cpp b/src/finish.cpp index d241adb..4a16d65 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -22,773 +22,796 @@ namespace hdrplus { - cv::Mat convert16bit2_8bit_(cv::Mat ans){ - if(ans.type()==CV_16UC3){ - cv::MatIterator_ it, end; - for( it = ans.begin(), end = ans.end(); it != end; ++it) - { - // std::cout< it, end; - for( it = ans.begin(), end = ans.end(); it != end; ++it) - { - // std::cout< it, end; - for( it = ans.begin(), end = ans.end(); it != end; ++it) - { - // std::cout<1){ - x = 1; - } - } - - x*=USHRT_MAX; - - return (uint16_t)x; - } - - uint16_t uGammaDecompress_1pix(float x, float threshold,float gainMin,float gainMax,float exponent){ - // Normalize pixel val - x/=65535.0; - // check the val against the threshold - if(x<=threshold){ - x = x/gainMin; - }else{ - x = pow((x+gainMax-1)/gainMax,exponent); - } - // clip - if(x<0){ - x=0; - }else{ - if(x>1){ - x = 1; - } - } - x*=65535; - - return (uint16_t)x; - } - - cv::Mat uGammaCompress_(cv::Mat m,float threshold,float gainMin,float gainMax,float exponent){ - if(m.type()==CV_16UC3){ - cv::MatIterator_ it, end; - for( it = m.begin(), end = m.end(); it != end; ++it) - { - (*it)[0] =uGammaCompress_1pix((*it)[0],threshold,gainMin,gainMax,exponent); - (*it)[1] =uGammaCompress_1pix((*it)[1],threshold,gainMin,gainMax,exponent); - (*it)[2] =uGammaCompress_1pix((*it)[2],threshold,gainMin,gainMax,exponent); - } - }else if(m.type()==CV_16UC1){ - uint16_t* ptr = (uint16_t*)m.data; - int end = m.rows*m.cols; - for(int i=0;i it, end; - for( it = m.begin(), end = m.end(); it != end; ++it) - { - (*it)[0] =uGammaDecompress_1pix((*it)[0],threshold,gainMin,gainMax,exponent); - (*it)[1] =uGammaDecompress_1pix((*it)[1],threshold,gainMin,gainMax,exponent); - (*it)[2] =uGammaDecompress_1pix((*it)[2],threshold,gainMin,gainMax,exponent); - } - }else if(m.type()==CV_16UC1){ - uint16_t* ptr = (uint16_t*)m.data; - int end = m.rows*m.cols; - for(int i=0;i it, end; + for (it = ans.begin(), end = ans.end(); it != end; ++it) + { + // std::cout< it, end; + for (it = ans.begin(), end = ans.end(); it != end; ++it) + { + // std::cout< it, end; + for (it = ans.begin(), end = ans.end(); it != end; ++it) + { + // std::cout< 1) { + x = 1; + } + } + + x *= USHRT_MAX; + + return (uint16_t)x; + } + + uint16_t uGammaDecompress_1pix(float x, float threshold, float gainMin, float gainMax, float exponent) { + // Normalize pixel val + x /= 65535.0; + // check the val against the threshold + if (x <= threshold) { + x = x / gainMin; + } + else { + x = pow((x + gainMax - 1) / gainMax, exponent); + } + // clip + if (x < 0) { + x = 0; + } + else { + if (x > 1) { + x = 1; + } + } + x *= 65535; + + return (uint16_t)x; + } + + cv::Mat uGammaCompress_(cv::Mat m, float threshold, float gainMin, float gainMax, float exponent) { + if (m.type() == CV_16UC3) { + cv::MatIterator_ it, end; + for (it = m.begin(), end = m.end(); it != end; ++it) + { + (*it)[0] = uGammaCompress_1pix((*it)[0], threshold, gainMin, gainMax, exponent); + (*it)[1] = uGammaCompress_1pix((*it)[1], threshold, gainMin, gainMax, exponent); + (*it)[2] = uGammaCompress_1pix((*it)[2], threshold, gainMin, gainMax, exponent); + } + } + else if (m.type() == CV_16UC1) { + uint16_t* ptr = (uint16_t*)m.data; + int end = m.rows*m.cols; + for (int i = 0; i < end; i++) { + *(ptr + i) = uGammaCompress_1pix(*(ptr + i), threshold, gainMin, gainMax, exponent); + } + + } + else { + std::cout << "Unsupported Data Type" << std::endl; + } + return m; + } + + cv::Mat uGammaDecompress_(cv::Mat m, float threshold, float gainMin, float gainMax, float exponent) { + if (m.type() == CV_16UC3) { + cv::MatIterator_ it, end; + for (it = m.begin(), end = m.end(); it != end; ++it) + { + (*it)[0] = uGammaDecompress_1pix((*it)[0], threshold, gainMin, gainMax, exponent); + (*it)[1] = uGammaDecompress_1pix((*it)[1], threshold, gainMin, gainMax, exponent); + (*it)[2] = uGammaDecompress_1pix((*it)[2], threshold, gainMin, gainMax, exponent); + } + } + else if (m.type() == CV_16UC1) { + uint16_t* ptr = (uint16_t*)m.data; + int end = m.rows*m.cols; + for (int i = 0; i < end; i++) { + *(ptr + i) = uGammaDecompress_1pix(*(ptr + i), threshold, gainMin, gainMax, exponent); + } + + } + else { + std::cout << "Unsupported Data Type" << std::endl; + } + + return m; + } + + cv::Mat gammasRGB(cv::Mat img, bool mode) { + if (mode) {// compress + return uGammaCompress_(img, 0.0031308, 12.92, 1.055, 1. / 2.4); + } + else { // decompress + return uGammaDecompress_(img, 0.04045, 12.92, 1.055, 2.4); + } + } + + void copy_mat_16U_2(uint16_t* ptr_A, cv::Mat B) { - // uint16_t* ptr_A = (uint16_t*)A.data; - uint16_t* ptr_B = (uint16_t*)B.data; - for(int r = 0; r < B.rows; r++) { - for(int c = 0; c < B.cols; c++) { - *(ptr_A+r*B.cols+c) = *(ptr_B+r*B.cols+c); - } - } - } - - cv::Mat mean_(cv::Mat img){ - // initialize processedImg - int H = img.rows; - int W = img.cols; - cv::Mat processedImg = cv::Mat(H,W,CV_16UC1); - uint16_t* ptr = (uint16_t*)processedImg.data; - - // traverse img - int idx = 0; - cv::MatIterator_ it, end; - for( it = img.begin(), end = img.end(); it != end; ++it) - { - uint32_t tmp = (*it)[0]+(*it)[1]+(*it)[2]; - uint16_t avg_val = tmp/3; - *(ptr+idx) = avg_val; - idx++; - } - - return processedImg; - } - - double getMean(cv::Mat img){ - uint16_t* ptr = (uint16_t*)img.data; - int max_idx = img.rows*img.cols*img.channels(); - double sum=0; - for(int i=0;iUSHRT_MAX){ - *(ptr+i) = USHRT_MAX; - }else{ - *(ptr+i)=(uint16_t)tmp; - } - } - return img; - } - - double getSaturated(cv::Mat img, double threshold){ - threshold *= USHRT_MAX; - double count=0; - uint16_t* ptr = (uint16_t*)img.data; - int max_idx = img.rows*img.cols*img.channels(); - for(int i=0;ithreshold){ - count++; - } - } - return count/(double)max_idx; - - } - - cv::Mat meanGain_(cv::Mat img,int gain){ - if(img.channels()!=3){ - std::cout<<"unsupport img type in meanGain_()"< it, end; - for( it = img.begin(), end = img.end(); it != end; ++it) - { - double sum = 0; - // R - double tmp = (*it)[0]*gain; - if(tmp<0) tmp=0; - if(tmp>USHRT_MAX) tmp = USHRT_MAX; - sum+=tmp; - - // G - tmp = (*it)[1]*gain; - if(tmp<0) tmp=0; - if(tmp>USHRT_MAX) tmp = USHRT_MAX; - sum+=tmp; - - // B - tmp = (*it)[2]*gain; - if(tmp<0) tmp=0; - if(tmp>USHRT_MAX) tmp = USHRT_MAX; - sum+=tmp; - - // put into processedImg - uint16_t avg_val = sum/3; - *(ptr+idx) = avg_val; - idx++; - } - return processedImg; - } - - } - - cv::Mat applyScaling_(cv::Mat mergedImage, cv::Mat shortGray, cv::Mat fusedGray){ - cv::Mat result = mergedImage.clone(); - uint16_t* ptr_shortg = (uint16_t*)shortGray.data; - uint16_t* ptr_fusedg = (uint16_t*)fusedGray.data; - int count = 0; - cv::MatIterator_ it, end; - for( it = result.begin(), end = result.end(); it != end; ++it) - { - double s = 1; - if(*(ptr_shortg+count)!=0){ - s = *(ptr_fusedg+count); - s/=*(ptr_shortg+count); - } - for(int c=0;cUSHRT_MAX){ - (*it)[c] = USHRT_MAX; - }else{ - (*it)[c] = tmp; - } - } - } - return result; - } - - void localToneMap(cv::Mat& mergedImage, Options options, cv::Mat& shortg, - cv::Mat& longg, cv::Mat& fusedg, int& gain){ - std::cout<<"HDR Tone Mapping..."<(mergedImage); //mean_(mergedImage); - std::cout<<"--- Compute grayscale image"< (1 - sSMean) / 2; // only works if burst underexposed - saturated = getSaturated(longSg,0.95); - if(options.verbose==4){ - - } - } - - }else{ - if(options.ltmGain>0){ - gain = options.ltmGain; - } - } - std::cout<<"--- Compute gain"< mergeMertens = cv::createMergeMertens(); - std::cout<<"--- Create Mertens"< src_expos; - src_expos.push_back(convert16bit2_8bit_(shortg.clone())); - src_expos.push_back(convert16bit2_8bit_(longg.clone())); - mergeMertens->process(src_expos, fusedg); - fusedg = fusedg*USHRT_MAX; - fusedg.convertTo(fusedg, CV_16UC1); - std::cout<<"--- Apply Mertens"<1){ - x = 1; - } - uint16_t result = x*USHRT_MAX; - return result; - } - - cv::Mat enhanceContrast(cv::Mat image, Options options){ - if(options.gtmContrast>=0 && options.gtmContrast<=1){ - uint16_t* ptr = (uint16_t*)image.data; - int end = image.rows*image.cols*image.channels(); - for(int idx = 0;idxUSHRT_MAX) r = USHRT_MAX; - *(ptr_r+idx) = (uint16_t)r; - } - return result; - } - - cv::Mat sharpenTriple(cv::Mat image, Tuning tuning, Options options){ - // sharpen the image using unsharp masking - std::vector amounts = tuning.sharpenAmount; - std::vector sigmas = tuning.sharpenSigma; - std::vector thresholds = tuning.sharpenThreshold; - // Compute all Gaussian blur - cv::Mat blur0,blur1,blur2; - cv::GaussianBlur(image,blur0,cv::Size(0,0),sigmas[0]); - cv::GaussianBlur(image,blur1,cv::Size(0,0),sigmas[1]); - cv::GaussianBlur(image,blur2,cv::Size(0,0),sigmas[2]); - std::cout<<" --- gaussian blur"< it, end; + for (it = img.begin(), end = img.end(); it != end; ++it) + { + uint32_t tmp = (*it)[0] + (*it)[1] + (*it)[2]; + uint16_t avg_val = tmp / 3; + *(ptr + idx) = avg_val; + idx++; + } + + return processedImg; + } + + double getMean(cv::Mat img) { + uint16_t* ptr = (uint16_t*)img.data; + int max_idx = img.rows*img.cols*img.channels(); + double sum = 0; + for (int i = 0; i < max_idx; i++) { + sum += *(ptr + i); + } + sum /= max_idx; + sum /= USHRT_MAX; + return sum; + } + + cv::Mat matMultiply_scalar(cv::Mat img, float gain) { + uint16_t* ptr = (uint16_t*)img.data; + int max_idx = img.rows*img.cols*img.channels(); + for (int i = 0; i < max_idx; i++) { + double tmp = *(ptr + i)*gain; + if (tmp < 0) { + *(ptr + i) = 0; + } + else if (tmp > USHRT_MAX) { + *(ptr + i) = USHRT_MAX; + } + else { + *(ptr + i) = (uint16_t)tmp; + } + } + return img; + } + + double getSaturated(cv::Mat img, double threshold) { + threshold *= USHRT_MAX; + double count = 0; + uint16_t* ptr = (uint16_t*)img.data; + int max_idx = img.rows*img.cols*img.channels(); + for (int i = 0; i < max_idx; i++) { + if (*(ptr + i) > threshold) { + count++; + } + } + return count / (double)max_idx; + + } + + cv::Mat meanGain_(cv::Mat img, int gain) { + if (img.channels() != 3) { + std::cout << "unsupport img type in meanGain_()" << std::endl; + return cv::Mat(); + } + else { // RGB img + int H = img.rows; + int W = img.cols; + cv::Mat processedImg = cv::Mat(H, W, CV_16UC1); + uint16_t* ptr = (uint16_t*)processedImg.data; + int idx = 0; + + cv::MatIterator_ it, end; + for (it = img.begin(), end = img.end(); it != end; ++it) + { + double sum = 0; + // R + double tmp = (*it)[0] * gain; + if (tmp < 0) tmp = 0; + if (tmp > USHRT_MAX) tmp = USHRT_MAX; + sum += tmp; + + // G + tmp = (*it)[1] * gain; + if (tmp < 0) tmp = 0; + if (tmp > USHRT_MAX) tmp = USHRT_MAX; + sum += tmp; + + // B + tmp = (*it)[2] * gain; + if (tmp < 0) tmp = 0; + if (tmp > USHRT_MAX) tmp = USHRT_MAX; + sum += tmp; + + // put into processedImg + uint16_t avg_val = sum / 3; + *(ptr + idx) = avg_val; + idx++; + } + return processedImg; + } + + } + + cv::Mat applyScaling_(cv::Mat mergedImage, cv::Mat shortGray, cv::Mat fusedGray) { + cv::Mat result = mergedImage.clone(); + uint16_t* ptr_shortg = (uint16_t*)shortGray.data; + uint16_t* ptr_fusedg = (uint16_t*)fusedGray.data; + int count = 0; + cv::MatIterator_ it, end; + for (it = result.begin(), end = result.end(); it != end; ++it) + { + double s = 1; + if (*(ptr_shortg + count) != 0) { + s = *(ptr_fusedg + count); + s /= *(ptr_shortg + count); + } + for (int c = 0; c < mergedImage.channels(); c++) { + double tmp = (*it)[c] * s; + if (tmp < 0) { + (*it)[c] = 0; + } + else if (tmp > USHRT_MAX) { + (*it)[c] = USHRT_MAX; + } + else { + (*it)[c] = tmp; + } + } + } + return result; + } + + void localToneMap(cv::Mat& mergedImage, Options options, cv::Mat& shortg, + cv::Mat& longg, cv::Mat& fusedg, int& gain) { + std::cout << "HDR Tone Mapping..." << std::endl; + // # Work with grayscale images + cv::Mat shortGray = rgb_2_gray(mergedImage); //mean_(mergedImage); + std::cout << "--- Compute grayscale image" << std::endl; + + // compute gain + gain = 0; + if (options.ltmGain == -1) { + double dsFactor = 25; + int down_height = round(shortGray.rows / dsFactor); + int down_width = round(shortGray.cols / dsFactor); + cv::Mat shortS; + cv::resize(shortGray, shortS, cv::Size(down_height, down_width), cv::INTER_LINEAR); + shortS = shortS.reshape(1, 1); + + bool bestGain = false; + double compression = 1.0; + double saturated = 0.0; + cv::Mat shortSg = gammasRGB(shortS.clone(), true); + double sSMean = getMean(shortSg); + + while ((compression < 1.9 && saturated < .95) || ((!bestGain) && (compression < 6) && (gain < 30) && (saturated < 0.33))) { + gain += 2; + cv::Mat longSg = gammasRGB(shortS.clone()*gain, true); + double lSMean = getMean(longSg); + compression = lSMean / sSMean; + bestGain = lSMean > (1 - sSMean) / 2; // only works if burst underexposed + saturated = getSaturated(longSg, 0.95); + if (options.verbose == 4) { + + } + } + + } + else { + if (options.ltmGain > 0) { + gain = options.ltmGain; + } + } + std::cout << "--- Compute gain" << std::endl; + // create a synthetic long exposure + cv::Mat longGray = meanGain_(mergedImage.clone(), gain); + std::cout << "--- Synthetic long expo" << std::endl; + // apply gamma correction to both + longg = gammasRGB(longGray.clone(), true); + shortg = gammasRGB(shortGray.clone(), true); + std::cout << "--- Apply Gamma correction" << std::endl; + // perform tone mapping by exposure fusion in grayscale + cv::Ptr mergeMertens = cv::createMergeMertens(); + std::cout << "--- Create Mertens" << std::endl; + // hack: cv2 mergeMertens expects inputs between 0 and 255 + // but the result is scaled between 0 and 1 (some values can actually be greater than 1!) + std::vector src_expos; + src_expos.push_back(convert16bit2_8bit_(shortg.clone())); + src_expos.push_back(convert16bit2_8bit_(longg.clone())); + mergeMertens->process(src_expos, fusedg); + fusedg = fusedg * USHRT_MAX; + fusedg.convertTo(fusedg, CV_16UC1); + std::cout << "--- Apply Mertens" << std::endl; + // undo gamma correction + cv::Mat fusedGray = gammasRGB(fusedg.clone(), false); + // cv::imwrite("fusedg_degamma.png", fusedGray); + std::cout << "--- Un-apply Gamma correction" << std::endl; + // scale each RGB channel of the short exposure accordingly + mergedImage = applyScaling_(mergedImage, shortGray, fusedGray); + std::cout << "--- Scale channels" << std::endl; + } + + uint16_t enhanceContrast_1pix(uint16_t pix_val, double gain) { + double x = pix_val; + x /= USHRT_MAX; + x = x - gain * sin(2 * M_PI*x); + if (x < 0) { + x = 0; + } + else if (x > 1) { + x = 1; + } + uint16_t result = x * USHRT_MAX; + return result; + } + + cv::Mat enhanceContrast(cv::Mat image, Options options) { + if (options.gtmContrast >= 0 && options.gtmContrast <= 1) { + uint16_t* ptr = (uint16_t*)image.data; + int end = image.rows*image.cols*image.channels(); + for (int idx = 0; idx < end; idx++) { + *(ptr + idx) = enhanceContrast_1pix(*(ptr + idx), options.gtmContrast); + } + } + else { + std::cout << "GTM ignored, expected a contrast enhancement ratio between 0 and 1" << std::endl; + } + return image; + } + + cv::Mat distL1_(cv::Mat X, cv::Mat Y) { + int end_x = X.rows*X.cols*X.channels(); + int end_y = Y.rows*Y.cols*Y.channels(); + cv::Mat result = cv::Mat(X.rows, X.cols, X.type()); + if (end_x == end_y) { + uint16_t* ptr_x = (uint16_t*)X.data; + uint16_t* ptr_y = (uint16_t*)Y.data; + uint16_t* ptr_r = (uint16_t*)result.data; + for (int i = 0; i < end_x; i++) { + if (*(ptr_x + i) < *(ptr_y + i)) { + *(ptr_r + i) = *(ptr_y + i) - *(ptr_x + i); + } + else { + *(ptr_r + i) = *(ptr_x + i) - *(ptr_y + i); + } + } + } + else { + std::cout << "Mat size not match. distL1_ failed!" << std::endl; + } + return result; + } + + cv::Mat sharpenTriple_(cv::Mat image, + cv::Mat blur0, cv::Mat low0, float th0, float k0, + cv::Mat blur1, cv::Mat low1, float th1, float k1, + cv::Mat blur2, cv::Mat low2, float th2, float k2) { + // create result mat + cv::Mat result = cv::Mat(image.rows, image.cols, image.type()); + // initialize iteraters + uint16_t* ptr_r = (uint16_t*)result.data; + uint16_t* ptr_img = (uint16_t*)image.data; + uint16_t* ptr_blur0 = (uint16_t*)blur0.data; + uint16_t* ptr_low0 = (uint16_t*)low0.data; + uint16_t* ptr_blur1 = (uint16_t*)blur1.data; + uint16_t* ptr_low1 = (uint16_t*)low1.data; + uint16_t* ptr_blur2 = (uint16_t*)blur2.data; + uint16_t* ptr_low2 = (uint16_t*)low2.data; + int n_channels = image.channels(); + int end = image.rows*image.cols*n_channels; + // traverse Image + for (int idx = 0; idx < end; idx++) { + double r, r0, r1, r2; + double x = *(ptr_img + idx); + double l0 = *(ptr_low0 + idx) / (double)USHRT_MAX; + double l1 = *(ptr_low1 + idx) / (double)USHRT_MAX; + double l2 = *(ptr_low2 + idx) / (double)USHRT_MAX; + double b0 = *(ptr_blur0 + idx); + double b1 = *(ptr_blur1 + idx); + double b2 = *(ptr_blur2 + idx); + r0 = l0 < th0 ? x : x + k0 * (x - b0); + r1 = l1 < th1 ? x : x + k1 * (x - b1); + r2 = l2 < th2 ? x : x + k2 * (x - b2); + r = (r0 + r1 + r2) / 3.0; + if (r < 0) r = 0; + if (r > USHRT_MAX) r = USHRT_MAX; + *(ptr_r + idx) = (uint16_t)r; + } + return result; + } + + cv::Mat sharpenTriple(cv::Mat image, Tuning tuning, Options options) { + // sharpen the image using unsharp masking + std::vector amounts = tuning.sharpenAmount; + std::vector sigmas = tuning.sharpenSigma; + std::vector thresholds = tuning.sharpenThreshold; + // Compute all Gaussian blur + cv::Mat blur0, blur1, blur2; + cv::GaussianBlur(image, blur0, cv::Size(0, 0), sigmas[0]); + cv::GaussianBlur(image, blur1, cv::Size(0, 0), sigmas[1]); + cv::GaussianBlur(image, blur2, cv::Size(0, 0), sigmas[2]); + std::cout << " --- gaussian blur" << std::endl; + // cv::imwrite("blur2.png", blur2); + // Compute all low contrast images + cv::Mat low0 = distL1_(blur0, image); + cv::Mat low1 = distL1_(blur1, image); + cv::Mat low2 = distL1_(blur2, image); + std::cout << " --- low contrast" << std::endl; + // cv::imwrite("low2.png", low2); + // Compute the triple sharpen + cv::Mat sharpImage = sharpenTriple_(image, + blur0, low0, thresholds[0], amounts[0], + blur1, low1, thresholds[1], amounts[1], + blur2, low2, thresholds[2], amounts[2]); + std::cout << " --- sharpen" << std::endl; + return sharpImage; + } + + void copy_mat_16U_3(uint16_t* ptr_A, cv::Mat B) { + // uint16_t* ptr_A = (uint16_t*)A.data; + uint16_t* ptr_B = (uint16_t*)B.data; + int H = B.rows; + int W = B.cols; + int end = H * W; + for (int i = 0; i < end; i++) { + *(ptr_A + i) = *(ptr_B + i); + } + } + + // void copy_mat_16U_3(uint16_t* ptr_A, cv::Mat B){ + // // uint16_t* ptr_A = (uint16_t*)A.data; + // uint16_t* ptr_B = (uint16_t*)B.data; + // for(int r = 0; r < B.rows; r++) { + // for(int c = 0; c < B.cols; c++) { + // *(ptr_A+r*B.cols+c) = *(ptr_B+r*B.cols+c); + // } + // } + // } + cv::Mat processMergedMat(cv::Mat mergedImg, int opencv_type) { + cv::Mat m; #if 0 - uint16_t* ptr = (uint16_t*)mergedImg.data; - for(int r = 0; r < mergedImg.rows; r++) { - std::vector dvals; - for(int c = 0; c < mergedImg.cols; c++) { - dvals.push_back(*(ptr+r*mergedImg.cols+c)); - } - cv::Mat mline(dvals, true); - cv::transpose(mline, mline); - m.push_back(mline); - } + uint16_t* ptr = (uint16_t*)mergedImg.data; + for (int r = 0; r < mergedImg.rows; r++) { + std::vector dvals; + for (int c = 0; c < mergedImg.cols; c++) { + dvals.push_back(*(ptr + r * mergedImg.cols + c)); + } + cv::Mat mline(dvals, true); + cv::transpose(mline, mline); + m.push_back(mline); + } #endif - int ch = CV_MAT_CN(opencv_type); - - m = mergedImg.clone(); - m = m.reshape(ch); - m.convertTo(m, opencv_type); - - return m; - } - - void show20_20(cv::Mat m){ - uint16_t* ptr = (uint16_t*)m.data; - for(int i=0;i<20;i++){ - for(int j=0;j<20;j++){ - std::cout<<*(ptr+i*m.cols+j)<<", "; - } - std::cout<refIdx = burst_images.reference_image_idx; - // this->burstPath = burstPath; - // std::cout<<"processMerged:"<refIdx = burst_images.reference_image_idx; + // this->burstPath = burstPath; + // std::cout<<"processMerged:"<mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); - // this->mergedBayer = processMergedMat(mergedB,CV_16UC1);//loadFromCSV("merged.csv", CV_16UC1); - // std::cout<<"processMerged:"<mergedBayer); - // this->mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); - // this->mergedBayer = processMergedMat(burst_images.merged_bayer_image, CV_16UC1); + this->mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); + // this->mergedBayer = processMergedMat(mergedB,CV_16UC1);//loadFromCSV("merged.csv", CV_16UC1); + // std::cout<<"processMerged:"<mergedBayer); + // this->mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); + // this->mergedBayer = processMergedMat(burst_images.merged_bayer_image, CV_16UC1); #else - // this->mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); - this->mergedBayer = processMergedMat(burst_images.merged_bayer_image, CV_16UC1); + // this->mergedBayer = loadFromCSV(DBG_OUTPUT_ROOT "merged.csv", CV_16UC1); + this->mergedBayer = processMergedMat(burst_images.merged_bayer_image, CV_16UC1); // std::cout<<"processMerged:"<mergedBayer); - // load_rawPathList(burstPath); + // std::cout<<"csv:"<mergedBayer); + // load_rawPathList(burstPath); // read in ref img - // bayer_image* ref = new bayer_image(rawPathList[refIdx]); - bayer_image* ref = new bayer_image(burst_images.bayer_images[burst_images.reference_image_idx]); - cv::Mat processedRefImage = postprocess(ref->libraw_processor,params.rawpyArgs); + // bayer_image* ref = new bayer_image(rawPathList[refIdx]); + bayer_image* ref = new bayer_image(burst_images.bayer_images[burst_images.reference_image_idx]); + cv::Mat processedRefImage = postprocess(ref->libraw_processor, params.rawpyArgs); - std::cout<<"size ref: "<refIdx]); - mergedImg->libraw_processor->imgdata.rawdata.raw_image = (uint16_t*)this->mergedBayer.data; - // copy_mat_16U_3(mergedImg->libraw_processor->imgdata.rawdata.raw_image,this->mergedBayer); - cv::Mat processedMerge = postprocess(mergedImg->libraw_processor,params.rawpyArgs); + // get the bayer_image of the merged image + // bayer_image* mergedImg = new bayer_image(rawPathList[refIdx]); + bayer_image* mergedImg = new bayer_image(burst_images.bayer_images[this->refIdx]); + mergedImg->libraw_processor->imgdata.rawdata.raw_image = (uint16_t*)this->mergedBayer.data; + // copy_mat_16U_3(mergedImg->libraw_processor->imgdata.rawdata.raw_image,this->mergedBayer); + cv::Mat processedMerge = postprocess(mergedImg->libraw_processor, params.rawpyArgs); -// write merged image + // write merged image #ifndef HDRPLUS_NO_DETAILED_OUTPUT - if(params.flags["writeMergedImage"]){ - std::cout<<"writing Merged img ..."<& libraw_ptr, cv::Mat B){ - uint16_t* ptr_A = (uint16_t*)libraw_ptr->imgdata.rawdata.raw_image; - uint16_t* ptr_B = (uint16_t*)B.data; - for(int r = 0; r < B.rows; r++) { - for(int c = 0; c < B.cols; c++) { - *(ptr_A+r*B.cols+c) = *(ptr_B+r*B.cols+c); - } - } - - } - - + // End of finishing + } + + void finish::copy_mat_16U(cv::Mat& A, cv::Mat B) { + uint16_t* ptr_A = (uint16_t*)A.data; + uint16_t* ptr_B = (uint16_t*)B.data; + for (int r = 0; r < A.rows; r++) { + for (int c = 0; c < A.cols; c++) { + *(ptr_A + r * A.cols + c) = *(ptr_B + r * B.cols + c); + } + } + } + + + + void finish::copy_rawImg2libraw(std::shared_ptr& libraw_ptr, cv::Mat B) { + uint16_t* ptr_A = (uint16_t*)libraw_ptr->imgdata.rawdata.raw_image; + uint16_t* ptr_B = (uint16_t*)B.data; + for (int r = 0; r < B.rows; r++) { + for (int c = 0; c < B.cols; c++) { + *(ptr_A + r * B.cols + c) = *(ptr_B + r * B.cols + c); + } + } + + } + + } // namespace hdrplus diff --git a/src/hdrplus_pipeline.cpp b/src/hdrplus_pipeline.cpp index 11e6fd2..8fa6a25 100644 --- a/src/hdrplus_pipeline.cpp +++ b/src/hdrplus_pipeline.cpp @@ -17,122 +17,122 @@ namespace hdrplus { -void hdrplus_pipeline::run_pipeline( \ - const std::string& burst_path, \ - const std::string& reference_image_path ) -{ - // Create burst of images - burst burst_images( burst_path, reference_image_path ); - - std::vector>>> alignments; - - // Run align - align_module.process( burst_images, alignments ); - - // Run merging - merge_module.process( burst_images, alignments ); - - - // Run finishing - cv::Mat finalImg; - finish_module.process( burst_images, finalImg); -} - -bool hdrplus_pipeline::run_pipeline( \ - const std::vector& burst_paths, \ - int reference_image_index, cv::Mat& finalImg ) -{ - // Create burst of images - burst burst_images( burst_paths, reference_image_index ); - std::vector>>> alignments; + void hdrplus_pipeline::run_pipeline(\ + const std::string& burst_path, \ + const std::string& reference_image_path) + { + // Create burst of images + burst burst_images(burst_path, reference_image_path); + + std::vector>>> alignments; + + // Run align + align_module.process(burst_images, alignments); + + // Run merging + merge_module.process(burst_images, alignments); + + + // Run finishing + cv::Mat finalImg; + finish_module.process(burst_images, finalImg); + } + + bool hdrplus_pipeline::run_pipeline(\ + const std::vector& burst_paths, \ + int reference_image_index, cv::Mat& finalImg) + { + // Create burst of images + burst burst_images(burst_paths, reference_image_index); + std::vector>>> alignments; #ifdef __ANDROID__ - // ALOGI("Finish loading images"); + // ALOGI("Finish loading images"); #endif - // Run align - align_module.process( burst_images, alignments ); + // Run align + align_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish align"); + // ALOGI("Finish align"); #endif - // Run merging - merge_module.process( burst_images, alignments ); + // Run merging + merge_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish merging"); + // ALOGI("Finish merging"); #endif - // Run finishing - finish_module.process( burst_images, finalImg); + // Run finishing + finish_module.process(burst_images, finalImg); #ifdef __ANDROID__ - // ALOGI("Finish process"); + // ALOGI("Finish process"); #endif - return true; -} + return true; + } -bool hdrplus_pipeline::run_pipeline( \ - const std::vector >& burst_contents, \ - int reference_image_index, cv::Mat& finalImg ) -{ - // Create burst of images - burst burst_images( burst_contents, reference_image_index ); - std::vector>>> alignments; + bool hdrplus_pipeline::run_pipeline(\ + const std::vector >& burst_contents, \ + int reference_image_index, cv::Mat& finalImg) + { + // Create burst of images + burst burst_images(burst_contents, reference_image_index); + std::vector>>> alignments; #ifdef __ANDROID__ - // ALOGI("Finish loading images"); + // ALOGI("Finish loading images"); #endif - // Run align - align_module.process( burst_images, alignments ); + // Run align + align_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish align"); + // ALOGI("Finish align"); #endif - // Run merging - merge_module.process( burst_images, alignments ); + // Run merging + merge_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish merging"); + // ALOGI("Finish merging"); #endif - // Run finishing - finish_module.process( burst_images, finalImg); + // Run finishing + finish_module.process(burst_images, finalImg); #ifdef __ANDROID__ - // ALOGI("Finish process"); + // ALOGI("Finish process"); #endif - return true; -} + return true; + } - bool hdrplus_pipeline::run_pipeline( \ - const std::vector >& burst_files, \ - int reference_image_index, cv::Mat& finalImg ) - { - // Create burst of images - burst burst_images( burst_files, reference_image_index ); - std::vector>>> alignments; + bool hdrplus_pipeline::run_pipeline(\ + const std::vector >& burst_files, \ + int reference_image_index, cv::Mat& finalImg) + { + // Create burst of images + burst burst_images(burst_files, reference_image_index); + std::vector>>> alignments; #ifdef __ANDROID__ - // ALOGI("Finish loading images"); + // ALOGI("Finish loading images"); #endif - // Run align - align_module.process( burst_images, alignments ); + // Run align + align_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish align"); + // ALOGI("Finish align"); #endif - // Run merging - merge_module.process( burst_images, alignments ); + // Run merging + merge_module.process(burst_images, alignments); #ifdef __ANDROID__ - // ALOGI("Finish merging"); + // ALOGI("Finish merging"); #endif - // Run finishing - finish_module.process( burst_images, finalImg); + // Run finishing + finish_module.process(burst_images, finalImg); #ifdef __ANDROID__ - // ALOGI("Finish process"); + // ALOGI("Finish process"); #endif - return true; - } + return true; + } } // namespace hdrplus diff --git a/src/merge.cpp b/src/merge.cpp index 9781c63..5f86776 100644 --- a/src/merge.cpp +++ b/src/merge.cpp @@ -8,345 +8,345 @@ namespace hdrplus { - void merge::process(hdrplus::burst& burst_images, \ - std::vector>>>& alignments) - { - // 4.1 Noise Parameters and RMS - // Noise parameters calculated from baseline ISO noise parameters - double lambda_shot, lambda_read; - std::tie(lambda_shot, lambda_read) = burst_images.bayer_images[burst_images.reference_image_idx].get_noise_params(); - - // 4.2-4.4 Denoising and Merging - - // Get padded bayer image - cv::Mat reference_image = burst_images.bayer_images_pad[burst_images.reference_image_idx]; + void merge::process(hdrplus::burst& burst_images, \ + std::vector>>>& alignments) + { + // 4.1 Noise Parameters and RMS + // Noise parameters calculated from baseline ISO noise parameters + double lambda_shot, lambda_read; + std::tie(lambda_shot, lambda_read) = burst_images.bayer_images[burst_images.reference_image_idx].get_noise_params(); + + // 4.2-4.4 Denoising and Merging + + // Get padded bayer image + cv::Mat reference_image = burst_images.bayer_images_pad[burst_images.reference_image_idx]; #ifndef NDEBUG - // cv::imwrite("ref.jpg", reference_image); + // cv::imwrite("ref.jpg", reference_image); #endif - // Get raw channels - std::vector channels(4); - hdrplus::extract_rgb_from_bayer(reference_image, channels[0], channels[1], channels[2], channels[3]); - - std::vector processed_channels(4); - // For each channel, perform denoising and merge - for (int i = 0; i < 4; ++i) { - // Get channel mat - cv::Mat channel_i = channels[i]; - // cv::imwrite("ref" + std::to_string(i) + ".jpg", channel_i); - - //we should be getting the individual channel in the same place where we call the processChannel function with the reference channel in its arguments - //possibly we could add another argument in the processChannel function which is the channel_i for the alternate image. maybe using a loop to cover all the other images - - //create list of channel_i of alternate images: - std::vector alternate_channel_i_list; - for (int j = 0; j < burst_images.num_images; j++) { - if (j != burst_images.reference_image_idx) { - - //get alternate image - cv::Mat alt_image = burst_images.bayer_images_pad[j]; - std::vector alt_channels(4); - hdrplus::extract_rgb_from_bayer(alt_image, alt_channels[0], alt_channels[1], alt_channels[2], alt_channels[3]); - - alternate_channel_i_list.push_back(alt_channels[i]); - } - } - - // Apply merging on the channel - cv::Mat merged_channel = processChannel(burst_images, alignments, channel_i, alternate_channel_i_list, lambda_shot, lambda_read); - // cv::imwrite("merged" + std::to_string(i) + ".jpg", merged_channel); - - // Put channel raw data back to channels - merged_channel.convertTo(processed_channels[i], CV_16U); - } - - // Write all channels back to a bayer mat - cv::Mat merged(reference_image.rows, reference_image.cols, CV_16U); - int x, y; - for (y = 0; y < reference_image.rows; ++y){ - uint16_t* row = merged.ptr(y); - if (y % 2 == 0){ - uint16_t* i0 = processed_channels[0].ptr(y / 2); - uint16_t* i1 = processed_channels[1].ptr(y / 2); - - for (x = 0; x < reference_image.cols;){ - //R - row[x] = i0[x / 2]; - x++; - - //G1 - row[x] = i1[x / 2]; - x++; - } - } - else { - uint16_t* i2 = processed_channels[2].ptr(y / 2); - uint16_t* i3 = processed_channels[3].ptr(y / 2); - - for(x = 0; x < reference_image.cols;){ - //G2 - row[x] = i2[x / 2]; - x++; - - //B - row[x] = i3[x / 2]; - x++; - } - } - } - - // Remove padding - std::vector padding = burst_images.padding_info_bayer; - cv::Range horizontal = cv::Range(padding[2], reference_image.cols - padding[3]); - cv::Range vertical = cv::Range(padding[0], reference_image.rows - padding[1]); - burst_images.merged_bayer_image = merged(vertical, horizontal); - // cv::imwrite("merged.jpg", burst_images.merged_bayer_image); - } - - void merge::getReferenceTiles(cv::Mat reference_image, std::vector& reference_tiles) { - for (int y = 0; y < reference_image.rows - offset; y += offset) { - for (int x = 0; x < reference_image.cols - offset; x += offset) { - cv::Mat tile = reference_image(cv::Rect(x, y, TILE_SIZE, TILE_SIZE)); - reference_tiles.push_back(tile); - } - } - } - - cv::Mat merge::mergeTiles(std::vector tiles, int num_rows, int num_cols) { - // 1. get all four subsets: original (evenly split), horizontal overlapped, - // vertical overlapped, 2D overlapped - std::vector> tiles_original; - std::vector row; - for (int y = 0; y < num_rows / offset - 1; y += 2) { - row.clear(); - for (int x = 0; x < num_cols / offset - 1; x += 2) { - row.push_back(tiles[y * (num_cols / offset - 1) + x]); - } - tiles_original.push_back(row); - } - - std::vector> tiles_horizontal; - // std::vector row; - for (int y = 0; y < num_rows / offset - 1; y += 2) { - row.clear(); - for (int x = 1; x < num_cols / offset - 1; x += 2) { - row.push_back(tiles[y * (num_cols / offset - 1) + x]); - } - tiles_horizontal.push_back(row); - } - - std::vector> tiles_vertical; - // std::vector row; - for (int y = 1; y < num_rows / offset - 1; y += 2) { - row.clear(); - for (int x = 0; x < num_cols / offset - 1; x += 2) { - row.push_back(tiles[y * (num_cols / offset - 1) + x]); - } - tiles_vertical.push_back(row); - } - - std::vector> tiles_2d; - // std::vector row; - for (int y = 1; y < num_rows / offset - 1; y += 2) { - row.clear(); - for (int x = 1; x < num_cols / offset - 1; x += 2) { - row.push_back(tiles[y * (num_cols / offset - 1) + x]); - } - tiles_2d.push_back(row); - } - - // 2. Concatenate the four subsets - cv::Mat img_original = cat2Dtiles(tiles_original); - cv::Mat img_horizontal = cat2Dtiles(tiles_horizontal); - cv::Mat img_vertical = cat2Dtiles(tiles_vertical); - cv::Mat img_2d = cat2Dtiles(tiles_2d); - - // 3. Add the four subsets together - img_original(cv::Rect(offset, 0, num_cols - TILE_SIZE, num_rows)) += img_horizontal; - img_original(cv::Rect(0, offset, num_cols, num_rows - TILE_SIZE)) += img_vertical; - img_original(cv::Rect(offset, offset, num_cols - TILE_SIZE, num_rows - TILE_SIZE)) += img_2d; - - return img_original; - } - - cv::Mat merge::processChannel(hdrplus::burst& burst_images, \ - std::vector>>>& alignments, \ - cv::Mat channel_image, \ - std::vector alternate_channel_i_list,\ - float lambda_shot, \ - float lambda_read) { - // Get tiles of the reference image + // Get raw channels + std::vector channels(4); + hdrplus::extract_rgb_from_bayer(reference_image, channels[0], channels[1], channels[2], channels[3]); + + std::vector processed_channels(4); + // For each channel, perform denoising and merge + for (int i = 0; i < 4; ++i) { + // Get channel mat + cv::Mat channel_i = channels[i]; + // cv::imwrite("ref" + std::to_string(i) + ".jpg", channel_i); + + //we should be getting the individual channel in the same place where we call the processChannel function with the reference channel in its arguments + //possibly we could add another argument in the processChannel function which is the channel_i for the alternate image. maybe using a loop to cover all the other images + + //create list of channel_i of alternate images: + std::vector alternate_channel_i_list; + for (int j = 0; j < burst_images.num_images; j++) { + if (j != burst_images.reference_image_idx) { + + //get alternate image + cv::Mat alt_image = burst_images.bayer_images_pad[j]; + std::vector alt_channels(4); + hdrplus::extract_rgb_from_bayer(alt_image, alt_channels[0], alt_channels[1], alt_channels[2], alt_channels[3]); + + alternate_channel_i_list.push_back(alt_channels[i]); + } + } + + // Apply merging on the channel + cv::Mat merged_channel = processChannel(burst_images, alignments, channel_i, alternate_channel_i_list, lambda_shot, lambda_read); + // cv::imwrite("merged" + std::to_string(i) + ".jpg", merged_channel); + + // Put channel raw data back to channels + merged_channel.convertTo(processed_channels[i], CV_16U); + } + + // Write all channels back to a bayer mat + cv::Mat merged(reference_image.rows, reference_image.cols, CV_16U); + int x, y; + for (y = 0; y < reference_image.rows; ++y) { + uint16_t* row = merged.ptr(y); + if (y % 2 == 0) { + uint16_t* i0 = processed_channels[0].ptr(y / 2); + uint16_t* i1 = processed_channels[1].ptr(y / 2); + + for (x = 0; x < reference_image.cols;) { + //R + row[x] = i0[x / 2]; + x++; + + //G1 + row[x] = i1[x / 2]; + x++; + } + } + else { + uint16_t* i2 = processed_channels[2].ptr(y / 2); + uint16_t* i3 = processed_channels[3].ptr(y / 2); + + for (x = 0; x < reference_image.cols;) { + //G2 + row[x] = i2[x / 2]; + x++; + + //B + row[x] = i3[x / 2]; + x++; + } + } + } + + // Remove padding + std::vector padding = burst_images.padding_info_bayer; + cv::Range horizontal = cv::Range(padding[2], reference_image.cols - padding[3]); + cv::Range vertical = cv::Range(padding[0], reference_image.rows - padding[1]); + burst_images.merged_bayer_image = merged(vertical, horizontal); + // cv::imwrite("merged.jpg", burst_images.merged_bayer_image); + } + + void merge::getReferenceTiles(cv::Mat reference_image, std::vector& reference_tiles) { + for (int y = 0; y < reference_image.rows - offset; y += offset) { + for (int x = 0; x < reference_image.cols - offset; x += offset) { + cv::Mat tile = reference_image(cv::Rect(x, y, TILE_SIZE, TILE_SIZE)); + reference_tiles.push_back(tile); + } + } + } + + cv::Mat merge::mergeTiles(std::vector tiles, int num_rows, int num_cols) { + // 1. get all four subsets: original (evenly split), horizontal overlapped, + // vertical overlapped, 2D overlapped + std::vector> tiles_original; + std::vector row; + for (int y = 0; y < num_rows / offset - 1; y += 2) { + row.clear(); + for (int x = 0; x < num_cols / offset - 1; x += 2) { + row.push_back(tiles[y * (num_cols / offset - 1) + x]); + } + tiles_original.push_back(row); + } + + std::vector> tiles_horizontal; + // std::vector row; + for (int y = 0; y < num_rows / offset - 1; y += 2) { + row.clear(); + for (int x = 1; x < num_cols / offset - 1; x += 2) { + row.push_back(tiles[y * (num_cols / offset - 1) + x]); + } + tiles_horizontal.push_back(row); + } + + std::vector> tiles_vertical; + // std::vector row; + for (int y = 1; y < num_rows / offset - 1; y += 2) { + row.clear(); + for (int x = 0; x < num_cols / offset - 1; x += 2) { + row.push_back(tiles[y * (num_cols / offset - 1) + x]); + } + tiles_vertical.push_back(row); + } + + std::vector> tiles_2d; + // std::vector row; + for (int y = 1; y < num_rows / offset - 1; y += 2) { + row.clear(); + for (int x = 1; x < num_cols / offset - 1; x += 2) { + row.push_back(tiles[y * (num_cols / offset - 1) + x]); + } + tiles_2d.push_back(row); + } + + // 2. Concatenate the four subsets + cv::Mat img_original = cat2Dtiles(tiles_original); + cv::Mat img_horizontal = cat2Dtiles(tiles_horizontal); + cv::Mat img_vertical = cat2Dtiles(tiles_vertical); + cv::Mat img_2d = cat2Dtiles(tiles_2d); + + // 3. Add the four subsets together + img_original(cv::Rect(offset, 0, num_cols - TILE_SIZE, num_rows)) += img_horizontal; + img_original(cv::Rect(0, offset, num_cols, num_rows - TILE_SIZE)) += img_vertical; + img_original(cv::Rect(offset, offset, num_cols - TILE_SIZE, num_rows - TILE_SIZE)) += img_2d; + + return img_original; + } + + cv::Mat merge::processChannel(hdrplus::burst& burst_images, \ + std::vector>>>& alignments, \ + cv::Mat channel_image, \ + std::vector alternate_channel_i_list, \ + float lambda_shot, \ + float lambda_read) { + // Get tiles of the reference image std::vector reference_tiles; getReferenceTiles(channel_image, reference_tiles); - // Get noise variance (sigma**2 = lambda_shot * tileRMS + lambda_read) - std::vector noise_variance = getNoiseVariance(reference_tiles, lambda_shot, lambda_read); - - // Apply FFT on reference tiles (spatial to frequency) - std::vector reference_tiles_DFT; - for (auto ref_tile : reference_tiles) { - cv::Mat ref_tile_DFT; - ref_tile.convertTo(ref_tile_DFT, CV_32F); - cv::dft(ref_tile_DFT, ref_tile_DFT, cv::DFT_COMPLEX_OUTPUT); - reference_tiles_DFT.push_back(ref_tile_DFT); - } - - // Acquire alternate tiles and apply FFT on them as well - std::vector> alt_tiles_list(reference_tiles.size()); - int num_tiles_row = alternate_channel_i_list[0].rows / offset - 1; - int num_tiles_col = alternate_channel_i_list[0].cols / offset - 1; - std::vector alt_tiles; - for (int y = 0; y < num_tiles_row; ++y) { - for (int x = 0; x < num_tiles_col; ++x) { - alt_tiles.clear(); - // Get reference tile location - int top_left_y = y * offset; - int top_left_x = x * offset; - - for (int i = 0; i < alternate_channel_i_list.size(); ++i) { - // Get alignment displacement - int displacement_y, displacement_x; - std::tie(displacement_y, displacement_x) = alignments[i + 1][y][x]; - // Get tile - cv::Mat alt_tile = alternate_channel_i_list[i](cv::Rect(top_left_x + displacement_x, top_left_y + displacement_y, TILE_SIZE, TILE_SIZE)); - // Apply FFT - cv::Mat alt_tile_DFT; - alt_tile.convertTo(alt_tile_DFT, CV_32F); - cv::dft(alt_tile_DFT, alt_tile_DFT, cv::DFT_COMPLEX_OUTPUT); - alt_tiles.push_back(alt_tile_DFT); - } - alt_tiles_list[y * num_tiles_col + x] = alt_tiles; - } - } - - // 4.2 Temporal Denoising - + // Get noise variance (sigma**2 = lambda_shot * tileRMS + lambda_read) + std::vector noise_variance = getNoiseVariance(reference_tiles, lambda_shot, lambda_read); + + // Apply FFT on reference tiles (spatial to frequency) + std::vector reference_tiles_DFT; + for (auto ref_tile : reference_tiles) { + cv::Mat ref_tile_DFT; + ref_tile.convertTo(ref_tile_DFT, CV_32F); + cv::dft(ref_tile_DFT, ref_tile_DFT, cv::DFT_COMPLEX_OUTPUT); + reference_tiles_DFT.push_back(ref_tile_DFT); + } + + // Acquire alternate tiles and apply FFT on them as well + std::vector> alt_tiles_list(reference_tiles.size()); + int num_tiles_row = alternate_channel_i_list[0].rows / offset - 1; + int num_tiles_col = alternate_channel_i_list[0].cols / offset - 1; + std::vector alt_tiles; + for (int y = 0; y < num_tiles_row; ++y) { + for (int x = 0; x < num_tiles_col; ++x) { + alt_tiles.clear(); + // Get reference tile location + int top_left_y = y * offset; + int top_left_x = x * offset; + + for (int i = 0; i < alternate_channel_i_list.size(); ++i) { + // Get alignment displacement + int displacement_y, displacement_x; + std::tie(displacement_y, displacement_x) = alignments[i + 1][y][x]; + // Get tile + cv::Mat alt_tile = alternate_channel_i_list[i](cv::Rect(top_left_x + displacement_x, top_left_y + displacement_y, TILE_SIZE, TILE_SIZE)); + // Apply FFT + cv::Mat alt_tile_DFT; + alt_tile.convertTo(alt_tile_DFT, CV_32F); + cv::dft(alt_tile_DFT, alt_tile_DFT, cv::DFT_COMPLEX_OUTPUT); + alt_tiles.push_back(alt_tile_DFT); + } + alt_tiles_list[y * num_tiles_col + x] = alt_tiles; + } + } + + // 4.2 Temporal Denoising + std::vector reference_tiles_DFTNew; - temporal_denoise(reference_tiles_DFT, alt_tiles_list, noise_variance, TEMPORAL_FACTOR, reference_tiles_DFTNew); + temporal_denoise(reference_tiles_DFT, alt_tiles_list, noise_variance, TEMPORAL_FACTOR, reference_tiles_DFTNew); reference_tiles_DFT.swap(reference_tiles_DFTNew); - // 4.3 Spatial Denoising + // 4.3 Spatial Denoising reference_tiles_DFTNew.clear(); - spatial_denoise(reference_tiles_DFT, alternate_channel_i_list.size(), noise_variance, SPATIAL_FACTOR, reference_tiles_DFTNew); + spatial_denoise(reference_tiles_DFT, alternate_channel_i_list.size(), noise_variance, SPATIAL_FACTOR, reference_tiles_DFTNew); reference_tiles_DFT.swap(reference_tiles_DFTNew); { std::vector empty; reference_tiles_DFTNew.swap(empty); } - //now reference tiles are temporally and spatially denoised - - // Apply IFFT on reference tiles (frequency to spatial) - std::vector denoised_tiles; - for (auto dft_tile : reference_tiles_DFT) { - cv::Mat denoised_tile; - cv::divide(dft_tile, TILE_SIZE * TILE_SIZE, dft_tile); - cv::dft(dft_tile, denoised_tile, cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT); - denoised_tiles.push_back(denoised_tile); - } + //now reference tiles are temporally and spatially denoised + + // Apply IFFT on reference tiles (frequency to spatial) + std::vector denoised_tiles; + for (auto dft_tile : reference_tiles_DFT) { + cv::Mat denoised_tile; + cv::divide(dft_tile, TILE_SIZE * TILE_SIZE, dft_tile); + cv::dft(dft_tile, denoised_tile, cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT); + denoised_tiles.push_back(denoised_tile); + } { reference_tiles.swap(denoised_tiles); std::vector empty; denoised_tiles.swap(empty); } - - - // 4.4 Cosine Window Merging - // Process tiles through 2D cosine window - std::vector windowed_tiles; - for (auto tile : reference_tiles) { - windowed_tiles.push_back(cosineWindow2D(tile)); - } - - // Merge tiles - return mergeTiles(windowed_tiles, channel_image.rows, channel_image.cols); - } - - void merge::temporal_denoise(std::vector tiles, std::vector> alt_tiles, std::vector noise_variance, float temporal_factor, std::vector& denoised) + + + // 4.4 Cosine Window Merging + // Process tiles through 2D cosine window + std::vector windowed_tiles; + for (auto tile : reference_tiles) { + windowed_tiles.push_back(cosineWindow2D(tile)); + } + + // Merge tiles + return mergeTiles(windowed_tiles, channel_image.rows, channel_image.cols); + } + + void merge::temporal_denoise(std::vector tiles, std::vector> alt_tiles, std::vector noise_variance, float temporal_factor, std::vector& denoised) { - // goal: temporially denoise using the weiner filter - // input: - // 1. array of 2D dft tiles of the reference image - // 2. array of 2D dft tiles of the aligned alternate image - // 3. estimated noise variance - // 4. temporal factor - // return: merged image patches dft - - // calculate noise scaling - double temporal_noise_scaling = (TILE_SIZE * TILE_SIZE * (2.0 / 16)) * TEMPORAL_FACTOR; - - // loop across tiles - for (int i = 0; i < tiles.size(); ++i) { - // sum of pairwise denoising - cv::Mat tile_sum = tiles[i].clone(); - double coeff = temporal_noise_scaling * noise_variance[i]; - - // Ref tile - cv::Mat tile = tiles[i]; - // Alt tiles - std::vector alt_tiles_i = alt_tiles[i]; - - for (int j = 0; j < alt_tiles_i.size(); ++j) { - // Alt tile - cv::Mat alt_tile = alt_tiles_i[j]; - // Tile difference - cv::Mat diff = tile - alt_tile; - - // Calculate absolute difference - cv::Mat complexMats[2]; - cv::split(diff, complexMats); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I)) - cv::magnitude(complexMats[0], complexMats[1], complexMats[0]); // planes[0] = magnitude - cv::Mat absolute_diff = complexMats[0].mul(complexMats[0]); - - // find shrinkage operator A - cv::Mat shrinkage; - cv::divide(absolute_diff, absolute_diff + coeff, shrinkage); - cv::merge(std::vector{shrinkage, shrinkage}, shrinkage); - - // Interpolation - tile_sum += alt_tile + diff.mul(shrinkage); - } - // Average by num of frames - cv::divide(tile_sum, alt_tiles_i.size() + 1, tile_sum); - denoised.push_back(tile_sum); - } - - - } - - void merge::spatial_denoise(std::vector tiles, int num_alts, std::vector noise_variance, float spatial_factor, std::vector& denoised) { - - double spatial_noise_scaling = (TILE_SIZE * TILE_SIZE * (1.0 / 16)) * spatial_factor; - - // Calculate |w| using ifftshift - cv::Mat row_distances = cv::Mat::zeros(1, TILE_SIZE, CV_32F); - for(int i = 0; i < TILE_SIZE; ++i) { - row_distances.at(i) = i - offset; - } - row_distances = cv::repeat(row_distances.t(), 1, TILE_SIZE); - cv::Mat col_distances = row_distances.t(); - cv::Mat distances; - cv::sqrt(row_distances.mul(row_distances) + col_distances.mul(col_distances), distances); - ifftshift(distances); - - // Loop through all tiles - for (int i = 0; i < tiles.size(); ++i) { - cv::Mat tile = tiles[i]; - float coeff = noise_variance[i] / (num_alts + 1) * spatial_noise_scaling; - - // Calculate absolute difference - cv::Mat complexMats[2]; - cv::split(tile, complexMats); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I)) - cv::magnitude(complexMats[0], complexMats[1], complexMats[0]); // planes[0] = magnitude - cv::Mat absolute_diff = complexMats[0].mul(complexMats[0]); - - // Division - cv::Mat scale; - cv::divide(absolute_diff, absolute_diff + distances * coeff, scale); - cv::merge(std::vector{scale, scale}, scale); - denoised.push_back(tile.mul(scale)); - } - } - + // goal: temporially denoise using the weiner filter + // input: + // 1. array of 2D dft tiles of the reference image + // 2. array of 2D dft tiles of the aligned alternate image + // 3. estimated noise variance + // 4. temporal factor + // return: merged image patches dft + + // calculate noise scaling + double temporal_noise_scaling = (TILE_SIZE * TILE_SIZE * (2.0 / 16)) * TEMPORAL_FACTOR; + + // loop across tiles + for (int i = 0; i < tiles.size(); ++i) { + // sum of pairwise denoising + cv::Mat tile_sum = tiles[i].clone(); + double coeff = temporal_noise_scaling * noise_variance[i]; + + // Ref tile + cv::Mat tile = tiles[i]; + // Alt tiles + std::vector alt_tiles_i = alt_tiles[i]; + + for (int j = 0; j < alt_tiles_i.size(); ++j) { + // Alt tile + cv::Mat alt_tile = alt_tiles_i[j]; + // Tile difference + cv::Mat diff = tile - alt_tile; + + // Calculate absolute difference + cv::Mat complexMats[2]; + cv::split(diff, complexMats); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I)) + cv::magnitude(complexMats[0], complexMats[1], complexMats[0]); // planes[0] = magnitude + cv::Mat absolute_diff = complexMats[0].mul(complexMats[0]); + + // find shrinkage operator A + cv::Mat shrinkage; + cv::divide(absolute_diff, absolute_diff + coeff, shrinkage); + cv::merge(std::vector{shrinkage, shrinkage}, shrinkage); + + // Interpolation + tile_sum += alt_tile + diff.mul(shrinkage); + } + // Average by num of frames + cv::divide(tile_sum, alt_tiles_i.size() + 1, tile_sum); + denoised.push_back(tile_sum); + } + + + } + + void merge::spatial_denoise(std::vector tiles, int num_alts, std::vector noise_variance, float spatial_factor, std::vector& denoised) { + + double spatial_noise_scaling = (TILE_SIZE * TILE_SIZE * (1.0 / 16)) * spatial_factor; + + // Calculate |w| using ifftshift + cv::Mat row_distances = cv::Mat::zeros(1, TILE_SIZE, CV_32F); + for (int i = 0; i < TILE_SIZE; ++i) { + row_distances.at(i) = i - offset; + } + row_distances = cv::repeat(row_distances.t(), 1, TILE_SIZE); + cv::Mat col_distances = row_distances.t(); + cv::Mat distances; + cv::sqrt(row_distances.mul(row_distances) + col_distances.mul(col_distances), distances); + ifftshift(distances); + + // Loop through all tiles + for (int i = 0; i < tiles.size(); ++i) { + cv::Mat tile = tiles[i]; + float coeff = noise_variance[i] / (num_alts + 1) * spatial_noise_scaling; + + // Calculate absolute difference + cv::Mat complexMats[2]; + cv::split(tile, complexMats); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I)) + cv::magnitude(complexMats[0], complexMats[1], complexMats[0]); // planes[0] = magnitude + cv::Mat absolute_diff = complexMats[0].mul(complexMats[0]); + + // Division + cv::Mat scale; + cv::divide(absolute_diff, absolute_diff + distances * coeff, scale); + cv::merge(std::vector{scale, scale}, scale); + denoised.push_back(tile.mul(scale)); + } + } + } // namespace hdrplus \ No newline at end of file