diff --git a/include/hdrplus/utility.h b/include/hdrplus/utility.h index ec13aca..8fca3c3 100644 --- a/include/hdrplus/utility.h +++ b/include/hdrplus/utility.h @@ -287,4 +287,59 @@ void extract_rgb_fmom_bayer( const cv::Mat& bayer_img, \ } } +/** + * @brief Extract RGB channel seprately from bayer image + * + * @tparam T data tyoe of bayer image. + * @return vector of RGB image. OpenCV internally maintain reference count. + * Thus this step won't create deep copy overhead. + * + * @example extract_rgb_fmom_bayer( bayer_img, rgb_vector_container ); + */ +template +void extract_rgb_fmom_bayer( const cv::Mat& bayer_img, \ + cv::Mat& red_img, cv::Mat& green_img1, cv::Mat& green_img2, cv::Mat& blue_img ) +{ + const T* bayer_img_ptr = (const T*)bayer_img.data; + int bayer_width = bayer_img.size().width; + int bayer_height = bayer_img.size().height; + int bayer_step = bayer_img.step1(); + + if ( bayer_width % 2 != 0 || bayer_height % 2 != 0 ) + { + throw std::runtime_error("Bayer image data size incorrect, must be multiplier of 2\n"); + } + + // RGB image is half the size of bayer image + int rgb_width = bayer_width / 2; + int rgb_height = bayer_height / 2; + red_img.create( rgb_height, rgb_width, bayer_img.type() ); + green_img1.create( rgb_height, rgb_width, bayer_img.type() ); + green_img2.create( rgb_height, rgb_width, bayer_img.type() ); + blue_img.create( rgb_height, rgb_width, bayer_img.type() ); + int rgb_step = red_img.step1(); + + T* r_img_ptr = (T*)red_img.data; + T* g1_img_ptr = (T*)green_img1.data; + T* g2_img_ptr = (T*)green_img2.data; + T* b_img_ptr = (T*)blue_img.data; + + for ( int rgb_row_i = 0; rgb_row_i < rgb_height; rgb_row_i++ ) + { + int rgb_row_i_offset = rgb_row_i * rgb_step; + + // Every RGB row corresbonding to two Bayer image row + int bayer_row_i_offset1 = ( rgb_row_i * 2 + 0 ) * bayer_step; // For RG + int bayer_row_i_offset2 = ( rgb_row_i * 2 + 1 ) * bayer_step; // For GB + + for ( int rgb_col_j = 0; rgb_col_j < rgb_width; rgb_col_j++ ) + { + r_img_ptr[ rgb_row_i_offset + rgb_col_j ] = bayer_img_ptr[ bayer_row_i_offset1 + ( rgb_col_j * 2 + 0 ) ]; + g1_img_ptr[ rgb_row_i_offset + rgb_col_j ] = bayer_img_ptr[ bayer_row_i_offset1 + ( rgb_col_j * 2 + 1 ) ]; + g2_img_ptr[ rgb_row_i_offset + rgb_col_j ] = bayer_img_ptr[ bayer_row_i_offset2 + ( rgb_col_j * 2 + 0 ) ]; + b_img_ptr[ rgb_row_i_offset + rgb_col_j ] = bayer_img_ptr[ bayer_row_i_offset2 + ( rgb_col_j * 2 + 1 ) ]; + } + } +} + } // namespace hdrplus diff --git a/src/merge.cpp b/src/merge.cpp index 95145d2..27495db 100644 --- a/src/merge.cpp +++ b/src/merge.cpp @@ -391,298 +391,4 @@ namespace hdrplus } -std::pair merge::getNoiseParams( int ISO, \ - int white_level, \ - double black_level ) -{ - // 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]; - - for (int y = 0; y < reference_image.rows; ++y) { - for (int x = 0; x < reference_image.cols; ++x) { - if (y % 2 == 0) { - if (x % 2 == 0) { - channels[0].push_back(reference_image.at(y, x)); - } else { - channels[1].push_back(reference_image.at(y, x)); - } - else { - if (x % 2 == 0) { - channels[2].push_back(reference_image.at(y, x)); - } - else { - channels[3].push_back(reference_image.at(y, x)); - } - } - } - } - - // For each channel, perform denoising and merge - for (int i = 0; i < 4; ++i) { - // Get channel mat - cv::Mat channel_i(reference_image.rows / 2, reference_image.cols / 2, CV_16U, channels[i].data()); - // cv::imwrite("ref" + std::to_string(i) + ".jpg", channel_i); - - // Apply merging on the channel - - //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_img_channel = getChannels(alt_image); //get channel array from alternate image - cv::Mat alt_channel_i(alt_image.rows / 2, alt_image.cols / 2, CV_16U, alt_img_channel[i].data()); - - alternate_channel_i_list.push_back(alt_channel_i) - } - } - - ///// - - //cv::Mat merged_channel = processChannel(burst_images, alignments, channel_i, lambda_shot, lambda_read); - - 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 - channels[i] = merged_channel.reshape(1, merged_channel.total()); - } - - // Write all channels back to a bayer mat - std::vector merged_raw; - - for (int y = 0; y < reference_image.rows; ++y) { - for (int x = 0; x < reference_image.cols; ++x) { - if (y % 2 == 0) { - if (x % 2 == 0) { - merged_raw.push_back(channels[0][(y / 2) * (reference_image.cols / 2) + (x / 2)]); - } - else { - merged_raw.push_back(channels[1][(y / 2) * (reference_image.cols / 2) + (x / 2)]); - } - } - else { - if (x % 2 == 0) { - merged_raw.push_back(channels[2][(y / 2) * (reference_image.cols / 2) + (x / 2)]); - } - else { - merged_raw.push_back(channels[3][(y / 2) * (reference_image.cols / 2) + (x / 2)]); - } - } - } - } - - // Create merged mat - cv::Mat merged(reference_image.rows, reference_image.cols, CV_16U, merged_raw.data()); - // cv::imwrite("merged.jpg", merged); - - // 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); - } - - 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); - } - - // TODO: 4.2 Temporal Denoising - - std::vector temporal_denoise(std::vector reference_tiles, std::vector reference_tiles_DFT, std::vector noise_varaince) { - //goal: temporially denoise using the weiner filter - //input: - //1. array of 2D dft tiles of the reference image - //2. array of 2D dft tiles ocf the aligned alternate image - //3. estimated noise varaince - //4. temporal factor - //return: merged image patches dft - - - - //tile_size = TILE_SIZE; - - double temporal_factor = 8.0 //8 by default - - double temporal_noise_scaling = (pow(TILE_SIZE,2) * (1.0/16*2))*temporal_factor; - - //start calculating the merged image tiles fft - - - //get the tiles of the alternate image as a list - - std::vector> alternate_channel_i_tile_list; //list of alt channel tiles - std::vector> alternate_tiles_DFT_list; //list of alt channel tiles - - for (auto alt_img_channel : alternate_channel_i_list) { - std::vector alt_img_channel_tile = getReferenceTiles(alt_img_channel); //get tiles from alt image - alternate_channel_i_tile_list.push_back(alt_img_channel_tile) - - std::vector alternate_tiles_DFT_list; - for (auto alt_tile : alt_img_channel_tile) { - 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); - alternate_tiles_DFT_list.push_back(alt_tile_DFT); - } - alternate_tiles_DFT_list.push_back(alternate_tiles_DFT); - } - - //get the dft of the alternate image - //std::vector alternate_tiles_DFT; - - - - - //std::vector tile_differences = reference_tiles_DFT - alternate_tiles_DFT_list; - - //find reference_tiles_DFT - alternate_tiles_DFT_list - std::vector> tile_difference_list; //list of tile differences - for (auto individual_alternate_tile_DFT : alternate_tiles_DFT_list) { - std::vector single_tile_difference = reference_tiles_DFT - individual_alternate_tile_DFT; - tile_difference_list.push_back(single_tile_difference); - } - - - // std::vector tile_sq_asolute_diff = tile_differences; //squared absolute difference is tile_differences.real**2 + tile_differnce.imag**2; //also tile_dist - - std::vector tile_sq_asolute_diff = tile_differences; //squared absolute difference is tile_differences.real**2 + tile_differnce.imag**2; //also tile_dist - - //get the real and imaginary components - /* - std::vector> absolute_difference_list; - for (auto individual_difference : tile_difference_list) { - for (int i =0; i < individual_difference.rows; i++ ) { - std::complex* row_ptr = tile_sq_asolute_diff.ptr>(i); - for (int j = 0; j< individual_difference.cols*individual_difference.channels(); j++) { - row_ptr = math.pow(individual_difference.at>(i,j).real(),2)+math.pow(individual_difference.at>(i,j).imag(),2); //.real and .imag - } - } - - //std::vector single_tile_difference = individual_difference.at>(0,0).real(); //.real and .imag - absolute_difference_list.push_back(single_tile_difference); - } - */ - - //find the squared absolute difference across all the tiles - - - std::vector A = tile_sq_asolute_diff/(tile_sq_asolute_diff+noise_variance) - - std::vector merged_image_tiles_fft = alternate_tiles_DFT_list + A * tile_differences; - - return merged_image_tiles_fft - - } - - std::vector spatial_denoise(std::vector reference_tiles, std::vector reference_tiles_DFT, std::vector noise_varaince) { - - double spatial_factor = 1; //to be added - double spatial_noise_scaling = (pow(TILE_SIZE,2) * (1.0/16*2))*spatial_factor; - - //calculate the spatial denoising - spatial_tile_dist = reference_tiles.real**2 + reference_tiles.imag**2; - std::vector WienerCoeff = denoised_tiles*spatial_noise_scaling*noise_variance; - - merged_channel_tiles_spatial = reference_tiles*spatial_tile_dist/(spatial_tile_dist+WienerCoeff) - - } - - } // namespace hdrplus \ No newline at end of file