diff --git a/include/hdrplus/burst.h b/include/hdrplus/burst.h index def6b78..c4c9dda 100644 --- a/include/hdrplus/burst.h +++ b/include/hdrplus/burst.h @@ -33,6 +33,9 @@ class burst // number of image (including reference) in burst int num_images; + + // Bayer image after merging, stored as cv::Mat + cv::Mat merged_bayer_image; }; } // namespace hdrplus diff --git a/include/hdrplus/merge.h b/include/hdrplus/merge.h index f980153..77319e7 100644 --- a/include/hdrplus/merge.h +++ b/include/hdrplus/merge.h @@ -2,14 +2,23 @@ #include #include // all opencv header +#include #include "hdrplus/burst.h" +#define TILE_SIZE 16 +#define TEMPORAL_FACTOR 75 +#define SPATIAL_FACTOR 0.1 + namespace hdrplus { class merge { public: + int offset = TILE_SIZE / 2; + float baseline_lambda_shot = 3.24 * pow( 10, -4 ); + float baseline_lambda_read = 4.3 * pow( 10, -6 ); + merge() = default; ~merge() = default; @@ -21,8 +30,155 @@ class merge * Outer most vector is per alternative image. * Inner most two vector is for horizontal & vertical tiles */ - void process( const hdrplus::burst& burst_images, \ + void process( hdrplus::burst& burst_images, \ std::vector>>>& alignments); + + + /* + std::vector get_other_tiles(); //return the other tile list T_1 to T_n + + std::vector vector_math(string operation, reference_tile, other_tile_list); //for loop allowing operations across single element and list + + std::vector scalar_vector_math(string operation, scalar num, std::vector tile_list); //for loop allowing operations across single element and list + + std::vector average_vector(std::vector tile_list); //take average of vector elements + + */ + + private: + float tileRMS(cv::Mat tile) { + cv::Mat squared; + cv::multiply(tile, tile, squared); + return sqrt(cv::mean(squared)[0]); + } + + std::vector getNoiseVariance(std::vector tiles, float lambda_shot, float lambda_read) { + std::vector noise_variance; + for (auto tile : tiles) { + noise_variance.push_back(lambda_shot * tileRMS(tile) + lambda_read); + } + return noise_variance; + } + + cv::Mat cosineWindow1D(cv::Mat input, int window_size = TILE_SIZE) { + cv::Mat output = input.clone(); + for (int i = 0; i < input.cols; ++i) { + output.at(0, i) = 1. / 2. - 1. / 2. * cos(2 * M_PI * (input.at(0, i) + 1 / 2.) / window_size); + } + return output; + } + + cv::Mat cosineWindow2D(cv::Mat tile) { + int window_size = tile.rows; // Assuming square tile + cv::Mat output_tile = tile.clone(); + + cv::Mat window = cv::Mat::zeros(1, window_size, CV_32F); + for(int i = 0; i < window_size; ++i) { + window.at(i) = i; + } + + cv::Mat window_x = cosineWindow1D(window, window_size); + window_x = cv::repeat(window_x, window_size, 1); + cv::Mat window_2d = window_x.mul(window_x.t()); + + cv::Mat window_applied; + cv::multiply(tile, window_2d, window_applied, 1, CV_32F); + return window_applied; + } + + cv::Mat cat2Dtiles(std::vector> tiles) { + std::vector rows; + for (auto row_tiles : tiles) { + cv::Mat row; + cv::hconcat(row_tiles, row); + rows.push_back(row); + } + cv::Mat img; + cv::vconcat(rows, img); + return img; + } + + void circshift(cv::Mat &out, const cv::Point &delta) + { + cv::Size sz = out.size(); + + // error checking + assert(sz.height > 0 && sz.width > 0); + + // no need to shift + if ((sz.height == 1 && sz.width == 1) || (delta.x == 0 && delta.y == 0)) + return; + + // delta transform + int x = delta.x; + int y = delta.y; + if (x > 0) x = x % sz.width; + if (y > 0) y = y % sz.height; + if (x < 0) x = x % sz.width + sz.width; + if (y < 0) y = y % sz.height + sz.height; + + + // in case of multiple dimensions + std::vector planes; + split(out, planes); + + for (size_t i = 0; i < planes.size(); i++) + { + // vertical + cv::Mat tmp0, tmp1, tmp2, tmp3; + cv::Mat q0(planes[i], cv::Rect(0, 0, sz.width, sz.height - y)); + cv::Mat q1(planes[i], cv::Rect(0, sz.height - y, sz.width, y)); + q0.copyTo(tmp0); + q1.copyTo(tmp1); + tmp0.copyTo(planes[i](cv::Rect(0, y, sz.width, sz.height - y))); + tmp1.copyTo(planes[i](cv::Rect(0, 0, sz.width, y))); + + // horizontal + cv::Mat q2(planes[i], cv::Rect(0, 0, sz.width - x, sz.height)); + cv::Mat q3(planes[i], cv::Rect(sz.width - x, 0, x, sz.height)); + q2.copyTo(tmp2); + q3.copyTo(tmp3); + tmp2.copyTo(planes[i](cv::Rect(x, 0, sz.width - x, sz.height))); + tmp3.copyTo(planes[i](cv::Rect(0, 0, x, sz.height))); + } + + cv::merge(planes, out); + } + + void fftshift(cv::Mat &out) + { + cv::Size sz = out.size(); + cv::Point pt(0, 0); + pt.x = (int) floor(sz.width / 2.0); + pt.y = (int) floor(sz.height / 2.0); + circshift(out, pt); + } + + void ifftshift(cv::Mat &out) + { + cv::Size sz = out.size(); + cv::Point pt(0, 0); + pt.x = (int) ceil(sz.width / 2.0); + pt.y = (int) ceil(sz.height / 2.0); + circshift(out, pt); + } + + std::vector getReferenceTiles(cv::Mat reference_image); + + cv::Mat mergeTiles(std::vector tiles, int rows, int cols); + + cv::Mat processChannel( hdrplus::burst& burst_images, \ + std::vector>>>& alignments, \ + cv::Mat channel_image, \ + std::vector alternate_channel_i_list,\ + float lambda_shot, \ + float lambda_read); + + //temporal denoise + std::vector temporal_denoise(std::vector tiles, std::vector> alt_tiles, std::vector noise_variance, float temporal_factor); + std::vector spatial_denoise(std::vector tiles, int num_alts, std::vector noise_variance, float spatial_factor); + + }; } // namespace hdrplus diff --git a/include/hdrplus/utility.h b/include/hdrplus/utility.h index 102f2ff..c192ab9 100644 --- a/include/hdrplus/utility.h +++ b/include/hdrplus/utility.h @@ -232,5 +232,4 @@ void print_img( const cv::Mat& img, int img_height = -1, int img_width = -1 ) printf("\n"); } - } // namespace hdrplus diff --git a/src/hdrplus_pipeline.cpp b/src/hdrplus_pipeline.cpp index 4508e1e..c2c0f68 100644 --- a/src/hdrplus_pipeline.cpp +++ b/src/hdrplus_pipeline.cpp @@ -27,7 +27,7 @@ void hdrplus_pipeline::run_pipeline( \ merge_module.process( burst_images, alignments ); // Run finishing - finish_module.process( burst_path, mergedBayer, refIdx); + finish_module.process( burst_path, burst_images.merged_bayer_image, burst_images.reference_image_idx); } } // namespace hdrplus diff --git a/src/merge.cpp b/src/merge.cpp index bb7b30b..0fa1e25 100644 --- a/src/merge.cpp +++ b/src/merge.cpp @@ -1,15 +1,333 @@ #include // all opencv header #include +#include #include "hdrplus/merge.h" #include "hdrplus/burst.h" +#include "hdrplus/utility.h" namespace hdrplus { -void merge::process( const hdrplus::burst& burst_images, \ - std::vector>>>& alignments) -{ + 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]; + cv::imwrite("ref.jpg", reference_image); + + // Get raw channels + std::vector channels(4); + hdrplus::extract_rgb_fmom_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_fmom_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); + } + + std::vector 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); + } + } + return reference_tiles; + } + + 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; + for (int y = 0; y < num_rows / offset - 1; y += 2) { + std::vector row; + 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; + for (int y = 0; y < num_rows / offset - 1; y += 2) { + std::vector row; + 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; + for (int y = 1; y < num_rows / offset - 1; y += 2) { + std::vector row; + 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; + for (int y = 1; y < num_rows / offset - 1; y += 2) { + std::vector row; + 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); + + // 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_SCALE | 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; + for (int y = 0; y < num_tiles_row; ++y) { + for (int x = 0; x < num_tiles_col; ++x) { + std::vector alt_tiles; + // 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_SCALE | 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 + reference_tiles_DFT = temporal_denoise(reference_tiles_DFT, alt_tiles_list, noise_variance, TEMPORAL_FACTOR); + + // 4.3 Spatial Denoising + reference_tiles_DFT = spatial_denoise(reference_tiles_DFT, alternate_channel_i_list.size(), noise_variance, SPATIAL_FACTOR); + //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::dft(dft_tile, denoised_tile, cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT); + denoised_tile.convertTo(denoised_tile, CV_16U); + denoised_tiles.push_back(denoised_tile); + } + reference_tiles = denoised_tiles; + + // 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); + } + + std::vector merge::temporal_denoise(std::vector tiles, std::vector> alt_tiles, std::vector noise_variance, float temporal_factor) { + // 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 = ((2.0 / 16)) * TEMPORAL_FACTOR; + + // loop across tiles + std::vector denoised; + for (int i = 0; i < tiles.size(); ++i) { + // sum of pairwise denoising + cv::Mat tile_sum = cv::Mat::zeros(TILE_SIZE, TILE_SIZE, CV_32FC2); + 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); + } + + return denoised; + } + + std::vector merge::spatial_denoise(std::vector tiles, int num_alts, std::vector noise_variance, float spatial_factor) { + + double spatial_noise_scaling = ((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); + + std::vector denoised; + // 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)); + } + return denoised; + } -} -} // namespace hdrplus +} // namespace hdrplus \ No newline at end of file