# -*- coding: utf-8; -*- """ Copyright (c) 2018 Rolf Hempel, rolf6419@gmx.de This file is part of the PlanetarySystemStacker tool (PSS). https://github.com/Rolf-Hempel/PlanetarySystemStacker PSS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with PSS. If not, see . """ from datetime import datetime from os import unlink from sys import stdout from time import time, sleep from cv2 import CV_32F, Laplacian, VideoWriter_fourcc, VideoWriter, FONT_HERSHEY_SIMPLEX, LINE_AA, \ putText, GaussianBlur, cvtColor, COLOR_BGR2HSV, COLOR_HSV2BGR, BORDER_DEFAULT, meanStdDev,\ resize, matchTemplate, minMaxLoc, TM_CCORR_NORMED, bilateralFilter, INTER_CUBIC from numpy import abs as np_abs from numpy import diff, average, hypot, sqrt, unravel_index, argmax, zeros, arange, array, matmul, \ empty, argmin, stack, sin, uint8, float32, uint16, full from math import exp from numpy import min as np_min from numpy.fft import fft2, ifft2 from numpy.linalg import solve from scipy.ndimage import sobel from exceptions import DivideByZeroError, ArgumentError, Error class Miscellaneous(object): """ This class provides static methods for various auxiliary purposes. """ @staticmethod def quality_measure(frame): """ Measure the amount of structure in a rectangular frame in both coordinate directions and return the minimum value. :param frame: 2D image :return: Measure for amount of local structure in the image (scalar) """ # Compute for each point the local gradient in both coordinate directions. dx = diff(frame) dy = diff(frame, axis=0) # Compute the sharpness per coordinate direction as the 1-norm of point values. sharpness_x = average(np_abs(dx)) sharpness_y = average(np_abs(dy)) # Return the sharpness in the direction where it is minimal. return min(sharpness_x, sharpness_y) @staticmethod def quality_measure_threshold(frame, black_threshold=40.): """ This is an alternative method for computing the amount of local structure. Here the summation only takes into account points where the luminosity exceeds a certain threshold. :param frame: 2D image :param black_threshold: Threshold for points to be considered :return: """ sum_horizontal = abs((frame[:, 2:] - frame[:, :-2])[frame[:, 1:-1] > black_threshold]).sum() sum_vertical = abs((frame[2:, :] - frame[:-2, :])[frame[1:-1, :] > black_threshold]).sum() return min(sum_horizontal, sum_vertical) @staticmethod def quality_measure_threshold_weighted(frame, stride=2, black_threshold=40., min_fraction=0.7): """ This is an alternative method for computing the amount of local structure. Here the summation only takes into account points where the luminosity exceeds a certain threshold. Additionally, if not too many points are discarded, the reduced sampling size is compensated. :param frame: 2D image :param stride: Stride for gradient computation. For blurry images increase this value. :param black_threshold: Threshold for points to be considered. :param min_fraction: Minimum fraction of points to pass the threshold. :return: """ frame_size = frame.shape[0] * frame.shape[1] stride_2 = 2 * stride # Compute a mask for all pixels which are bright enough (to avoid background noise). mask = frame > black_threshold mask_fraction = mask.sum() / frame_size # If most pixels are bright enough, compensate for different pixel counts. if mask_fraction > min_fraction: sum_horizontal = abs((frame[:, stride_2:] - frame[:, :-stride_2])[mask[:, stride:-stride]]).sum() / mask_fraction sum_vertical = abs((frame[stride_2:, :] - frame[:-stride_2, :])[mask[stride:-stride, :]]).sum() / mask_fraction # If many pixels are too dim, penalize this patch by not compensating for pixel count. else: sum_horizontal = abs((frame[:, stride_2:] - frame[:, :-stride_2])[mask[:, stride:-stride]]).sum() sum_vertical = abs((frame[stride_2:, :] - frame[:-stride_2, :])[mask[stride:-stride, :]]).sum() return min(sum_horizontal, sum_vertical) @staticmethod def local_contrast_laplace(frame, stride): """ Measure the amount of structure in a rectangular frame in both coordinate directions and return the minimum value. The discrete variance of the Laplacian is computed. :param frame: 2D image :param stride: Factor for down-sampling :return: Measure for amount of local structure in the image (scalar) """ # sharpness = sum(laplace(frame[::stride, ::stride])**2) return meanStdDev(Laplacian(frame[::stride, ::stride], CV_32F))[1][0][0] @staticmethod def local_contrast_sobel(frame, stride): """ Compute a measure for local contrast in an image using the Sobel method. :param frame: 2D image :param stride: Factor for down-sampling :return: Overall sharpness measure (scalar) """ frame_int32 = frame[::stride, ::stride].astype('int32') dx = sobel(frame_int32, 0) # vertical derivative dy = sobel(frame_int32, 1) # horizontal derivative mag = hypot(dx, dy) # magnitude return mag.sum(axis=0) @staticmethod def local_contrast(frame, stride): """ Compute a measure for local contrast in an image based on local gradients. Down-sample the image by factor "stride" before measuring the contrast. :param frame: 2D image :param stride: Factor for down-sampling :return: Overall measure for local contrast (scalar) """ frame_strided = frame[::stride, ::stride] # Remove a row or column, respectively, to make the dx and dy arrays of the same shape. dx = diff(frame_strided)[1:, :] # remove the first row dy = diff(frame_strided, axis=0)[:, 1:] # remove the first column dnorm = hypot(dx, dy) return average(dnorm) @staticmethod def translation(frame_0, frame_1, shape): """ Compute the translation vector of frame_1 relative to frame_0 using phase correlation. :param frame_0: Reference frame :param frame_1: Frame shifted slightly relative to frame_0 :param shape: Shape of frames :return: [shift in y, shift in x] of frame_1 relative to frame_0. More precisely: frame_1 must be shifted by this amount to register with frame_0. """ # Compute the fast Fourier transforms of both frames. f0 = fft2(frame_0) f1 = fft2(frame_1) # Compute the phase correlation. The resulting image has a bright peak at the pixel location # which corresponds to the shift vector. ir = abs(ifft2((f0 * f1.conjugate()) / (abs(f0) * abs(f1)))) # Compute the pixel coordinates of the image maximum. ty, tx = unravel_index(argmax(ir), shape) # Bring the shift values as close as possible to the coordinate origin. if ty > shape[0] // 2: ty -= shape[0] if tx > shape[1] // 2: tx -= shape[1] return [ty, tx] @staticmethod def multilevel_correlation(reference_box_first_phase, frame_mono_blurred, blurr_strength_first_phase, reference_box_second_phase, y_low, y_high, x_low, x_high, search_width, weight_matrix_first_phase=None, subpixel_solve=False): """ Determine the local warp shift at an alignment point using a multi-level approach based on normalized cross correlation. The first level uses a pixel grid which is coarser by a factor of two in both coordinate directions. The second level uses the original pixel grid. An additional noise reduction is applied on the first level as chosen by the corresponding blurring parameter. The global search width is split between the two phases: A fixed search width of 4 is used in the second phase (only for local corrections). Therefore, a width of (search_width-4) in each coordinate direction remains for the first phase. In both phases it is determined if the optimum is attained on the border of the search area. If so, the correlation is regarded as unsuccessful (because the real optimum may be outside of the search area). :param reference_box_first_phase: Image box with stride 2 around alignment point in the locally sharpest frame. A Gaussian filter with strength "blurr_strength_first_phase" has been applied. :param frame_mono_blurred: Given frame (stride 1) for which the local shift at the alignment point is to be computed. This is the Gaussian blurred version of the monochrome frame as computed in class "frames". :param blurr_strength_first_phase: Additional Gaussian blur strength to be applied to images in the first phase. :param reference_box_second_phase: Image box with stride 1 around alignment point in the locally sharpest frame. :param y_low: Lower y coordinate limit of box in given frame, taking into account the global shift and the different sizes of the mean frame and the original frames. :param y_high: Upper y coordinate limit. :param x_low: Lower x coordinate limit. :param x_high: Upper x coordinate limit. :param search_width: Maximum distance in y and x from origin of the search area. (See the comment in the explanatory text above.) :param weight_matrix_first_phase: This parameter (if not None) may give a weighting array by which the cross correlation results are multiplied before the maximum value is determined. The size of this 2D array in each coordinate direction is that of the reference_box_first_phase plus two times the first phase search width (see below). :param subpixel_solve: If True, in the second phase the optimum is computed with sub-pixel accuracy (i.e. the returned shifts are not integer). If False, shifts are computed as integer values. :return: (shift_y_local_first_phase, shift_x_local_first_phase, success_first_phase, shift_y_local_second_phase, shift_x_local_second_phase, success_second_phase) with: shift_y_local_first_phase: Local y warp shift determined in first phase (expressed in terms of the original pixel grid). shift_x_local_first_phase: Local x warp shift determined in first phase. success_first_phase: "True" if the optimum was attained in the interior of the search domain. "False" otherwise. shift_y_local_second_phase: Local y warp shift determined in second phase. shift_x_local_second_phase: Local x warp shift determined in second phase. success_second_phase: "True" if the optimum was attained in the interior of the search domain. "False" otherwise. """ # The optimization in the second phase is only meant for local fine grid adjustments. search_width_second_phase = 4 # In the first phase the largest part of the local warp is detected. Divide the search # width by two because the first phase uses a coarser pixel grid. search_width_first_phase = int((search_width - search_width_second_phase) / 2) # Define a window around the alignment point box. Extend the window in each coordinate # direction. This defines the search space for the template matching. Coarsen the grid # by a factor of two and apply an additional Gaussian blur. index_extension = search_width_first_phase * 2 frame_window_first_phase = GaussianBlur(frame_mono_blurred[ y_low - index_extension:y_high + index_extension:2, x_low - index_extension:x_high + index_extension:2], (blurr_strength_first_phase, blurr_strength_first_phase), 0) # Compute the normalized cross correlation. result = matchTemplate((frame_window_first_phase).astype(float32), reference_box_first_phase.astype(float32), TM_CCORR_NORMED) # Determine the position of the maximum correlation and compute the corresponding warp # shift. The factor of 2 transforms the shift to the fine pixel grid. If a non-trivial # weight matrix is specified, multiply the results before looking for the maximum position. if weight_matrix_first_phase is not None: minVal, maxVal, minLoc, maxLoc = minMaxLoc(result * weight_matrix_first_phase) else: minVal, maxVal, minLoc, maxLoc = minMaxLoc(result) shift_y_local_first_phase = (search_width_first_phase - maxLoc[1]) * 2 shift_x_local_first_phase = (search_width_first_phase - maxLoc[0]) * 2 # The first phase is regarded as successful if the maximum correlation was attained in the # interior of the search space. success_first_phase = abs(shift_y_local_first_phase) != index_extension and abs( shift_x_local_first_phase) != index_extension # If the first phase was successful, add a second phase with a local search on the finest # (original) pixel grid. if success_first_phase: # Define the search window for the second phase. Apply the shift found in the # first phase. y_lo = y_low - shift_y_local_first_phase - search_width_second_phase y_hi = y_high - shift_y_local_first_phase + search_width_second_phase x_lo = x_low - shift_x_local_first_phase - search_width_second_phase x_hi = x_high - shift_x_local_first_phase + search_width_second_phase # Cut out the frame window for the second phase (fine grid) correlation. frame_window_second_phase = frame_mono_blurred[y_lo:y_hi, x_lo:x_hi] # Perform the template matching on the fine grid, again using normalized cross # correlation. result = matchTemplate(frame_window_second_phase.astype(float32), reference_box_second_phase, TM_CCORR_NORMED) # Find the position of the local optimum and compute the corresponding shift values. minVal, maxVal, minLoc, maxLoc = minMaxLoc(result) shift_y_local_second_phase = search_width_second_phase - maxLoc[1] shift_x_local_second_phase = search_width_second_phase - maxLoc[0] # Again, the second phase is deemed successful only if the optimum was attained in the # interior of the search space. success_second_phase = abs( shift_y_local_second_phase) != search_width_second_phase and abs( shift_x_local_second_phase) != search_width_second_phase # If the second phase was not successful, set the corresponding shifts to zero. if not success_second_phase: shift_y_local_second_phase = shift_x_local_second_phase = 0 # The following code computes the sub-pixel shift correction to be used in drizzling. elif subpixel_solve: # Cut a 3x3 window around the optimum from the matching results. surroundings = result[maxLoc[1]-1:maxLoc[1]+2, maxLoc[0]-1:maxLoc[0]+2] try: # Compute the correction to the center position. Only if the found sub-pixel # correction is within the 3x3 box, it is trusted (and used). y_corr, x_corr = Miscellaneous.sub_pixel_solve(surroundings) if abs(y_corr) <= 1. and abs(x_corr) <= 1.: shift_y_local_second_phase -= y_corr shift_x_local_second_phase -= x_corr except: # print ("Subpixel solve not successful") pass # If the first phase was unsuccessful, drop the second phase and set all warp shifts to 0. else: success_second_phase = False shift_y_local_first_phase = shift_x_local_first_phase = shift_y_local_second_phase = \ shift_x_local_second_phase = 0 return shift_y_local_first_phase, shift_x_local_first_phase, success_first_phase, \ shift_y_local_second_phase, shift_x_local_second_phase, success_second_phase @staticmethod def search_local_match(reference_box, frame, y_low, y_high, x_low, x_high, search_width, sampling_stride, sub_pixel=True): """ Try shifts in y, x between the box around the alignment point in the mean frame and the corresponding box in the given frame. Start with shifts [0, 0] and move out in a circular fashion, until the radius "search_width" is reached. The global frame shift is accounted for beforehand already. :param reference_box: Image box around alignment point in mean frame. :param frame: Given frame for which the local shift at the alignment point is to be computed. :param y_low: Lower y coordinate limit of box in given frame, taking into account the global shift and the different sizes of the mean frame and the original frames. :param y_high: Upper y coordinate limit :param x_low: Lower x coordinate limit :param x_high: Upper x coordinate limit :param search_width: Maximum radius of the search spiral :param sampling_stride: Stride in both coordinate directions used in computing deviations :param sub_pixel: If True, compute local shifts with sub-pixel accuracy :return: ([shift_y, shift_x], [min_r]) with: shift_y, shift_x: shift values of minimum or [0, 0] if no optimum could be found. [dev_r]: list of minimum deviations for radius r=0 ... r_max, where r_max is the widest search radius tested. """ # Initialize the global optimum with an impossibly large value. deviation_min = 1.e30 dy_min = None dx_min = None # Initialize list of minimum deviations for each search radius and field of deviations. dev_r = [] deviations = zeros((2 * search_width + 1, 2 * search_width + 1)) # Start with shift [0, 0] and proceed in a circular pattern. for r in arange(search_width + 1): # Create an enumerator which produces shift values [dy, dx] in a circular pattern # with radius "r". circle_r = Miscellaneous.circle_around(0, 0, r) # Initialize the optimum for radius "r" to an impossibly large value, # and the corresponding shifts to None. deviation_min_r, dy_min_r, dx_min_r = 1.e30, None, None # Go through the circle with radius "r" and compute the difference (deviation) # between the shifted frame and the corresponding box in the mean frame. Find the # minimum "deviation_min_r" for radius "r". if sampling_stride != 1: for (dy, dx) in circle_r: deviation = abs( reference_box[::sampling_stride, ::sampling_stride] - frame[ y_low - dy:y_high - dy:sampling_stride, x_low - dx:x_high - dx:sampling_stride]).sum() # deviation = sqrt(square( # reference_box - frame[y_low - dy:y_high - dy, x_low - dx:x_high - dx]).sum()) if deviation < deviation_min_r: deviation_min_r, dy_min_r, dx_min_r = deviation, dy, dx deviations[dy + search_width, dx + search_width] = deviation else: for (dy, dx) in circle_r: deviation = abs( reference_box - frame[y_low - dy:y_high - dy, x_low - dx:x_high - dx]).sum() if deviation < deviation_min_r: deviation_min_r, dy_min_r, dx_min_r = deviation, dy, dx deviations[dy + search_width, dx + search_width] = deviation # Append the minimal deviation for radius r to list of minima. dev_r.append(deviation_min_r) # If for the current radius there is no improvement compared to the previous radius, # the optimum is reached. if deviation_min_r >= deviation_min: # For sub-pixel accuracy, find local minimum of fitting paraboloid. if sub_pixel: try: y_correction, x_correction = Miscellaneous.sub_pixel_solve(deviations[ dy_min + search_width - 1: dy_min + search_width + 2, dx_min + search_width - 1: dx_min + search_width + 2]) except DivideByZeroError as ex: print(ex.message) x_correction = y_correction = 0. # Add the sub-pixel correction to the local shift. If the correction is above # one pixel, something went wrong. if abs(y_correction) < 1. and abs(x_correction) < 1.: dy_min += y_correction dx_min += x_correction # else: # print ("y_correction: " + str(y_correction) + ", x_correction: " + # str(x_correction)) return [dy_min, dx_min], dev_r # Otherwise, update the current optimum and continue. else: deviation_min, dy_min, dx_min = deviation_min_r, dy_min_r, dx_min_r # If within the maximum search radius no optimum could be found, return [0, 0]. return [0, 0], dev_r # Define the matrix used in method "sub_pixel_solve" to solve the normal equations. It is # computed once and for all as "inv(a_transpose * a) * a_transpose". Using this matrix, the # solution of the normal equations is reduced to a simple matrix multiplication. sub_pixel_solve_matrix = [ [0.16666667, -0.33333333, 0.16666667, 0.16666667, -0.33333333, 0.16666667, 0.16666667, -0.33333333, 0.16666667], [0.16666667, 0.16666667, 0.16666667, -0.33333333, -0.33333333, -0.33333333, 0.16666667, 0.16666667, 0.16666667], [0.25, 0., -0.25, 0., 0., 0., -0.25, 0., 0.25], [-0.16666667, 0., 0.16666667, -0.16666667, 0., 0.16666667, -0.16666667, 0., 0.16666667], [-0.16666667, -0.16666667, -0.16666667, 0., 0., 0., 0.16666667, 0.16666667, 0.16666667], [-0.11111111, 0.22222222, -0.11111111, 0.22222222, 0.55555556, 0.22222222, -0.11111111, 0.22222222, -0.11111111]] @staticmethod def sub_pixel_solve(function_values): """ Compute the sub-pixel correction for method "search_local_match". :param function_values: Matching differences at (3 x 3) pixels around the minimum / maximum found :return: Corrections in y and x to the center position for local minimum / maximum """ # If the functions are not yet reduced to 1D, do it now. function_values_1d = function_values.reshape((9,)) # Solve for parameters of the fitting function: # f = a_f * x ** 2 + b_f * y ** 2 + c_f * x * y + d_f * x + e_f * y + g_f # using normal equations. The problem is reduced to a matrix multiplication with the fixed # system matrix defined above. a_f, b_f, c_f, d_f, e_f, g_f = matmul(Miscellaneous.sub_pixel_solve_matrix, function_values_1d) # print("\nSolve, coeffs: " + str((a_f, b_f, c_f, d_f, e_f, g_f))) # The corrected pixel values of the minimum result from setting the first derivatives of # the fitting funtion in y and x direction to zero, and solving for y and x. denominator_y = c_f ** 2 - 4. * a_f * b_f if abs(denominator_y) > 1.e-10 and abs(a_f) > 1.e-10: y_correction = (2. * a_f * e_f - c_f * d_f) / denominator_y x_correction = (- c_f * y_correction - d_f) / (2. * a_f) elif abs(denominator_y) > 1.e-10 and abs(c_f) > 1.e-10: y_correction = (2. * a_f * e_f - c_f * d_f) / denominator_y x_correction = (-2. * b_f * y_correction - e_f) / c_f else: raise DivideByZeroError("Sub-pixel shift cannot be computed, set to zero") return y_correction, x_correction @staticmethod def sub_pixel_solve_old(function_values): """ Compute the sub-pixel correction for method "search_local_match". :param function_values: Matching differences at (3 x 3) pixels around the minimum found :return: Corrections in y and x to the center position for local minimum """ # If the functions are not yet reduced to 1D, do it now. function_values_1d = function_values.reshape((9,)) # There are nine equations for six unknowns. Use normal equations to solve for optimum. a_transpose = array( [[1., 0., 1., 1., 0., 1., 1., 0., 1.], [1., 1., 1., 0., 0., 0., 1., 1., 1.], [1., -0., -1., -0., 0., 0., -1., 0., 1.], [-1., 0., 1., -1., 0., 1., -1., 0., 1.], [-1., -1., -1., 0., 0., 0., 1., 1., 1.], [1., 1., 1., 1., 1., 1., 1., 1., 1.]]) a_transpose_a = array( [[6., 4., 0., 0., 0., 6.], [4., 6., 0., 0., 0., 6.], [0., 0., 4., 0., 0., 0.], [0., 0., 0., 6., 0., 0.], [0., 0., 0., 0., 6., 0.], [6., 6., 0., 0., 0., 9.]]) # Right hand side is "a transposed times input vector". rhs = matmul(a_transpose, function_values_1d) # Solve for parameters of the fitting function # f = a_f * x ** 2 + b_f * y ** 2 + c_f * x * y + d_f * x + e_f * y + g_f a_f, b_f, c_f, d_f, e_f, g_f = solve(a_transpose_a, rhs) # The corrected pixel values of the minimum result from setting the first derivatives of # the fitting funtion in y and x direction to zero, and solving for y and x. denominator_y = c_f ** 2 - 4. * a_f * b_f if abs(denominator_y) > 1.e-10 and abs(a_f) > 1.e-10: y_correction = (2. * a_f * e_f - c_f * d_f) / denominator_y x_correction = (- c_f * y_correction - d_f) / (2. * a_f) elif abs(denominator_y) > 1.e-10 and abs(c_f) > 1.e-10: y_correction = (2. * a_f * e_f - c_f * d_f) / denominator_y x_correction = (-2. * b_f * y_correction - e_f) / c_f else: raise DivideByZeroError("Sub-pixel shift cannot be computed, set to zero") return y_correction, x_correction @staticmethod def search_local_match_init(reference_frame, y_low, y_high, x_low, x_high, search_width): """ As an alternative to the method "search_local_match" above there is a split version consisting of this initialization method and the execution method. For a given alignment point the initialization has to be called once only. The data structures allocated during initialization are re-used with each frame for which the shift against the reference frame is to be computed. The individual frame shift is computed with the execution method below. The split version can be better suited for execution on a graphic accelerator. The biggest part of the data does not change between frame executions, so it has to be copied to the graphics memory only once. The executions for individual frames only need little additional data, and they can be performed in parallel. The split version is not available with subpixel accuracy. :param reference_frame: Reference frame, against which the current frame shift is to be computed (note: The full frame, not the alignment box) :param y_low: Lower y coordinate limit of alignment box in reference frame :param y_high: Upper y coordinate limit of alignment box in reference frame :param x_low: Lower x coordinate limit of alignment box in reference frame :param x_high: Upper x coordinate limit of alignment box in reference frame :param search_width: Maximum radius of the search spiral :return: (reference_stack, displacements, radius_start) witn reference_stack: Stack of shifted reference frame alignment boxes displacements: List of [dy, dx] displacements for all stack elements radius_start: Start indices into "reference_stack" and "displacements" for each radius (from 0 to search_width) """ # Compute the maximal number of shifts to be tested. search_dim = (2*search_width+1)**2 # Allocate permanent data structures. window_height = y_high - y_low window_width = x_high - x_low reference_stack = empty([search_dim, window_height, window_width], dtype=reference_frame.dtype) displacements = [] radius_start = [0] # Create the reference_stack containing shifted windows into the reference frame. index = 0 for r in range(search_width + 1): # Create an enumerator which produces shift values [dy, dx] in a circular pattern # with radius "r". for (dy, dx) in Miscellaneous.circle_around(0, 0, r): reference_stack[index, :, :] = reference_frame[y_low + dy:y_high + dy, x_low + dx:x_high + dx] displacements.append([dy, dx]) index += 1 radius_start.append(index) return (reference_stack, displacements, radius_start) @staticmethod def search_local_match_execute(frame_window, reference_stack, displacements, radius_start): """ The "execute" method of the split version of "search_local_match". This method is executed for each frame, potentially in parallel for all frames. :param frame_window: Window into the current frame of the same shape as specified for the reference window in the init method. Be careful: the current frame and the reference frame have different shapes. :param reference_stack: Permanent object allocated by init routing (see above) :param displacements: Permanent object allocated by init routing (see above) :param radius_start: Permanent object allocated by init routing (see above) :return: ([shift_y, shift_x], [min_r]) with: shift_y, shift_x: shift values of minimum or [0, 0] if no optimum could be found. [dev_r]: list of minimum deviations for radius r=0 ... r_max, where r_max is the widest search radius tested. """ search_width_plus_1 = len(radius_start) # Initialize the global optimum with an impossibly large value. deviation_min = 1.e30 index_min = None # Initialize list of minimum deviations for each search radius and field of deviations. dev_r = [] # Compare frame_window with a stack of shifted reference frame windows stored for radius r. for r in arange(search_width_plus_1): temp_vec = abs( reference_stack[radius_start[r]:radius_start[r + 1], :, :] - frame_window).sum( axis=(1, 2)) deviation_min_r = np_min(temp_vec) index_min_r = argmin(temp_vec) + radius_start[r] # The same in loop notation: # # deviation_min_r, index_min_r = 1.e30, None # for index in arange(radius_start[r], radius_start[r+1]): # deviation = abs(reference_stack[index, :,:] - frame_window).sum() # if deviation < deviation_min_r: # deviation_min_r = deviation # index_min_r = index # Append the minimal deviation for radius r to list of minima. dev_r.append(deviation_min_r) if deviation_min_r >= deviation_min: return displacements[index_min], dev_r # Otherwise, update the current optimum and continue. else: deviation_min = deviation_min_r index_min = index_min_r # If within the maximum search radius no optimum could be found, return [0, 0]. return [0, 0], dev_r @staticmethod def search_local_match_gradient(reference_box, frame, y_low, y_high, x_low, x_high, search_width, sampling_stride, dev_table): """ Try shifts in y, x between the box around the alignment point in the mean frame and the corresponding box in the given frame. Start with shifts [0, 0] and move out in steps until a local optimum is reached. In each step try all positions with distance 1 in y and/or x from the optimum found in the previous step (steepest descent). The global frame shift is accounted for beforehand already. :param reference_box: Image box around alignment point in mean frame. :param frame: Given frame for which the local shift at the alignment point is to be computed. :param y_low: Lower y coordinate limit of box in given frame, taking into account the global shift and the different sizes of the mean frame and the original frames. :param y_high: Upper y coordinate limit :param x_low: Lower x coordinate limit :param x_high: Upper x coordinate limit :param search_width: Maximum distance in y and x from origin of the search area :param sampling_stride: Stride in both coordinate directions used in computing deviations :param dev_table: Scratch table to be used internally for storing intermediate results, size: [2*search_width, 2*search_width], dtype=float32. :return: ([shift_y, shift_x], [min_r]) with: shift_y, shift_x: shift values of minimum or [0, 0] if no optimum could be found. [dev_r]: list of minimum deviations for all steps until a local minimum is found. """ # Set up a table which keeps deviation values from earlier iteration steps. This way, # deviation evaluations can be avoided at coordinates which have been visited before. # Initialize deviations with an impossibly high value. dev_table[:, :] = 1.e30 # Initialize the global optimum with the value at dy=dx=0. if sampling_stride != 1: deviation_min = abs(reference_box[::sampling_stride, ::sampling_stride] - frame[ y_low:y_high:sampling_stride, x_low:x_high:sampling_stride]).sum() else: deviation_min = abs(reference_box - frame[y_low:y_high, x_low:x_high]).sum() dev_table[0, 0] = deviation_min dy_min = 0 dx_min = 0 counter_new = 0 counter_reused = 0 # Initialize list of minimum deviations for each search radius. dev_r = [deviation_min] # Start with shift [0, 0]. Stop when a circle with radius 1 around the current optimum # reaches beyond the search area. while max(abs(dy_min), abs(dx_min)) < search_width-1: # Create an enumerator which produces shift values [dy, dx] in a circular pattern # with radius 1 around the current optimum [dy_min, dx_min]. circle_1 = Miscellaneous.circle_around(dy_min, dx_min, 1) # Initialize the optimum for the new circle to an impossibly large value, # and the corresponding shifts to None. deviation_min_1, dy_min_1, dx_min_1 = 1.e30, None, None # Go through the circle with radius 1 and compute the difference (deviation) # between the shifted frame and the corresponding box in the mean frame. Find the # minimum "deviation_min_1". if sampling_stride != 1: for (dy, dx) in circle_1: deviation = dev_table[dy, dx] if deviation > 1.e29: counter_new += 1 deviation = abs(reference_box[::sampling_stride, ::sampling_stride] - frame[ y_low - dy:y_high - dy:sampling_stride, x_low - dx:x_high - dx:sampling_stride]).sum() dev_table[dy, dx] = deviation else: counter_reused += 1 if deviation < deviation_min_1: deviation_min_1, dy_min_1, dx_min_1 = deviation, dy, dx # deviation = abs(reference_box[::sampling_stride, ::sampling_stride] - frame[ # y_low - dy:y_high - dy:sampling_stride, # x_low - dx:x_high - dx:sampling_stride]).sum() # if deviation < deviation_min_1: # deviation_min_1, dy_min_1, dx_min_1 = deviation, dy, dx else: for (dy, dx) in circle_1: deviation = dev_table[dy, dx] if deviation > 1.e29: deviation = abs( reference_box - frame[y_low - dy:y_high - dy, x_low - dx:x_high - dx]).sum() dev_table[dy, dx] = deviation if deviation < deviation_min_1: deviation_min_1, dy_min_1, dx_min_1 = deviation, dy, dx # Append the minimal deviation found in this step to list of minima. dev_r.append(deviation_min_1) # If for the current center the match is better than for all neighboring points, a # local optimum is found. if deviation_min_1 >= deviation_min: # print ("new: " + str(counter_new) + ", reused: " + str(counter_reused)) return [dy_min, dx_min], dev_r # Otherwise, update the current optimum and continue. else: deviation_min, dy_min, dx_min = deviation_min_1, dy_min_1, dx_min_1 # If within the maximum search radius no optimum could be found, return [0, 0]. return [0, 0], dev_r @staticmethod def search_local_match_full(reference_box, frame, y_low, y_high, x_low, x_high, search_width, sampling_stride, dev_table): """ For all shifts with -search_width < shift < search_width in y, x between the box around the alignment point in the mean frame and the corresponding box in the given frame compute the deviation. Return the y and x offsets where the deviation is smallest. The global frame shift is accounted for beforehand already. :param reference_box: Image box around alignment point in mean frame. :param frame: Given frame for which the local shift at the alignment point is to be computed. :param y_low: Lower y coordinate limit of box in given frame, taking into account the global shift and the different sizes of the mean frame and the original frames. :param y_high: Upper y coordinate limit :param x_low: Lower x coordinate limit :param x_high: Upper x coordinate limit :param search_width: Maximum distance in y and x from origin of the search area :param sampling_stride: Stride in both coordinate directions used in computing deviations :param dev_table: Scratch table to be used internally for storing intermediate results, size: [2*search_width+1, 2*search_width+1], dtype=float32. :return: ([shift_y, shift_x], [min_r]) with: shift_y, shift_x: shift values of minimum or [0, 0] if no optimum could be found. [dev_r]: list of minimum deviations for all steps until a local minimum is found. """ # Set up a table which keeps deviation values from earlier iteration steps. This way, # deviation evaluations can be avoided at coordinates which have been visited before. # Initialize deviations with an impossibly high value. dev_table[:, :] = 1.e30 if sampling_stride != 1: for dy in range(-search_width, search_width + 1): for dx in range(-search_width, search_width + 1): dev_table[search_width + dy, search_width + dx] = abs( reference_box[::sampling_stride, ::sampling_stride] - frame[ y_low - dy:y_high - dy:sampling_stride, x_low - dx:x_high - dx:sampling_stride]).sum() else: for dy in range(-search_width, search_width + 1): for dx in range(-search_width, search_width + 1): dev_table[search_width + dy, search_width + dx] = abs( reference_box - frame[y_low - dy:y_high - dy, x_low - dx:x_high - dx]).sum() dy, dx = unravel_index(argmin(dev_table, axis=None), dev_table.shape) # Return the coordinate shifts for the minimum position and the value at that position. return [dy-search_width, dx-search_width], dev_table[dy, dx] @staticmethod def auto_rgb_align(input_image, max_shift, interpolation_factor=1, reduce_output=True, blur_strength=None): """ Align the three color channels of an RGB image automatically. For sub-pixel resolution the image can be interpolated before the shift is measured. Optionally, a Gaussian blur can be added to suppress noise. If no matching shift can be found within the range given by "max_shift", an Exception of type Error is thrown. :param input_image: Three-channel RGB image. :param max_shift: Maximal displacement between channels. :param interpolation_factor: Scaling factor (integer) for subpixel measurements. :param reduce_output: If True, the corrected image is reduced to the input resolution. If False, it stays at the interpolated resolution. In this case the output values for correction_red and correction_blue are in terms of the interpolated resolution as well. :param blur_strength: Optional blur strength, must be an uneven integer > 0. :return: (corrected_image, correction_red, correction_blue) with: corrected_image: The corrected image with the same datatype as the input image. Please note that the size may be reduced because of channel mismatch at the borders. correction_red: Tuple (shift_y, shift_x) with coordinate shifts in y and x applied to the red channel of input_image to produce the corrected_image. correction_blue: Tuple (shift_y, shift_x) with coordinate shifts in y and x applied to the blue channel of input_image to produce the corrected_image. """ # sleep(2.) # Immediately return for monochrome input. if len(input_image.shape) != 3: return input_image # If subpixel resolution is asked for, interpolate the input image first. if interpolation_factor != 1: input_interpolated = Miscellaneous.shift_colors(input_image, (0, 0), (0, 0), interpolate_input=interpolation_factor) else: input_interpolated = input_image # Measure the shifts of the red and blue channels, respectively, with respect to the green # channel. channel_green = 1 channel_red = 0 shift_red = Miscellaneous.measure_rgb_shift(input_interpolated, channel_red, channel_green, max_shift * interpolation_factor, blur_strength=blur_strength) channel_blue = 2 shift_blue = Miscellaneous.measure_rgb_shift(input_interpolated, channel_blue, channel_green, max_shift*interpolation_factor, blur_strength=blur_strength) # Reverse the shift measured in the input image. if reduce_output: factor = interpolation_factor else: factor = 1 return Miscellaneous.shift_colors(input_interpolated, (-shift_red[0], -shift_red[1]), (-shift_blue[0], -shift_blue[1]), reduce_output=factor), \ (-shift_red[0]/factor, -shift_red[1]/factor), \ (-shift_blue[0]/factor, -shift_blue[1]/factor) @staticmethod def shift_colors(input_image, shift_red, shift_blue, interpolate_input=1, reduce_output=1): """ Shift the red and blue channel of a color image in y and x direction, and leave the green channel unchanged. The shift is sub-pixel accurate if the input is interpolated. Optionally the resulting image is reduced in size by a given factor. :param input_image: Three-channel RGB image. :param shift_red: Tuple (shift_y, shift_x) with shifts in y and x direction for the red channel. After multiplication with the interpolate_input factor (if given) the shifts are rounded to the next integer value. This enables sub-pixel shifts at the original image scale. :param shift_blue: Tuple (shift_y, shift_x) with shifts in y and x direction for the blue channel. :param interpolate_input: If set, it must be an integer >= 1. The input image is interpolated in both directions by this factor before the shift is applied. :param reduce_output: If set, it must be an integer >= 1. The intermediate image is reduced in size by this factor. Please note that interpolate_input = reduce_output assures that the input and output images are the same size. :return: Three-channel RGB image with the color shifts applied. """ # sleep(1.) # Immediately return for monochrome input. if len(input_image.shape) != 3: return input_image # If all shifts are zero, nothing is to be done except interpolation / reduction. if not (shift_red[0] or shift_red[1] or shift_blue[0] or shift_blue[1]): # If the factors for interpolation and reduction are the same, nothing is to be done. if not (interpolate_input - reduce_output) or not reduce_output: return input_image # Simple resizing. reduce_output is not zero (see above)! else: scale_factor = float(interpolate_input) / float(reduce_output) return resize(input_image, (round(input_image.shape[1] * scale_factor), round(input_image.shape[0] * scale_factor)), interpolation=INTER_CUBIC) # If interpolation is requested, resize the input image and multiply the shift values. if interpolate_input != 1: dim_y = input_image.shape[0] * interpolate_input dim_x = input_image.shape[1] * interpolate_input interp_image = resize(input_image, (dim_x, dim_y), interpolation=INTER_CUBIC) s_red_y = round(interpolate_input * shift_red[0]) s_red_x = round(interpolate_input * shift_red[1]) s_blue_y = round(interpolate_input * shift_blue[0]) s_blue_x = round(interpolate_input * shift_blue[1]) else: dim_y = input_image.shape[0] dim_x = input_image.shape[1] interp_image = input_image s_red_y = round(shift_red[0]) s_red_x = round(shift_red[1]) s_blue_y = round(shift_blue[0]) s_blue_x = round(shift_blue[1]) # Remember the data type of the input image. The output image will be of the same type. type = input_image.dtype # In the following, for both the red and blue channels the index areas are computed for # which no shift is beyond the image dimensions. First treat the y coordinate. y_low_source_r = -s_red_y y_high_source_r = dim_y - s_red_y y_low_target_r = 0 # If the shift reaches beyond the frame, reduce the copy area. if y_low_source_r < 0: y_low_target_r = -y_low_source_r y_low_source_r = 0 if y_high_source_r > dim_y: y_high_source_r = dim_y y_high_target_r = y_low_target_r + y_high_source_r - y_low_source_r y_low_source_b = -s_blue_y y_high_source_b = dim_y - s_blue_y y_low_target_b = 0 # If the shift reaches beyond the frame, reduce the copy area. if y_low_source_b < 0: y_low_target_b = -y_low_source_b y_low_source_b = 0 if y_high_source_b > dim_y: y_high_source_b = dim_y y_high_target_b = y_low_target_b + y_high_source_b - y_low_source_b # The same for the x coordinate. x_low_source_r = -s_red_x x_high_source_r = dim_x - s_red_x x_low_target_r = 0 # If the shift reaches beyond the frame, reduce the copy area. if x_low_source_r < 0: x_low_target_r = -x_low_source_r x_low_source_r = 0 if x_high_source_r > dim_x: x_high_source_r = dim_x x_high_target_r = x_low_target_r + x_high_source_r - x_low_source_r x_low_source_b = -s_blue_x x_high_source_b = dim_x - s_blue_x x_low_target_b = 0 # If the shift reaches beyond the frame, reduce the copy area. if x_low_source_b < 0: x_low_target_b = -x_low_source_b x_low_source_b = 0 if x_high_source_b > dim_x: x_high_source_b = dim_x x_high_target_b = x_low_target_b + x_high_source_b - x_low_source_b # Now the coordinate bounds can be computed for which the output image has entries in all # three channels. The green channel is not shifted and can just be copied in place. min_y_g = max(y_low_target_r, y_low_target_b) max_y_g = min(y_high_target_r, y_high_target_b) min_x_g = max(x_low_target_r, x_low_target_b) max_x_g = min(x_high_target_r, x_high_target_b) # Allocate space for the output image (still not reduced in size), and copy the three color # channels with the appropriate shifts applied. output_image_interp = empty((max_y_g - min_y_g, max_x_g - min_x_g, 3), dtype=type) output_image_interp[:, :, 0] = interp_image[min_y_g - s_red_y:max_y_g - s_red_y, min_x_g - s_red_x:max_x_g - s_red_x, 0] output_image_interp[:, :, 1] = interp_image[min_y_g:max_y_g, min_x_g:max_x_g, 1] output_image_interp[:, :, 2] = interp_image[min_y_g - s_blue_y:max_y_g - s_blue_y, min_x_g - s_blue_x:max_x_g - s_blue_x, 2] # If output size reduction is specified, resize the intermediate image. if reduce_output != 1: output_image = resize(output_image_interp, (round(output_image_interp.shape[1] / reduce_output), round(output_image_interp.shape[0] / reduce_output)), interpolation=INTER_CUBIC) else: output_image = output_image_interp return output_image @staticmethod def measure_rgb_shift(image, channel_id, reference_id, max_shift, blur_strength=None): """ Measure the shift between two color channels of a three-color RGB image. Before measuring the shift with cross correlation, a Gaussian blur can be applied to both channels to reduce noise. :param image: Input image (3 channel RGB). :param channel_id: Index of the channel, the shift of which is to be measured against the reference channel. :param reference_id: Index of the reference channel, usually 1 for "green". :param max_shift: Maximal search space radius in pixels. :param blur_strength: Strength of the Gaussian blur to be applied first. :return: Tuple (shift_y, shift_x) with coordinate shifts in y and x. """ # Test for invalid channel ids. if not (0 <= channel_id <= 2): raise ArgumentError("Invalid color channel id: " + str(channel_id)) if not (0 <= reference_id <= 2): raise ArgumentError("Invalid reference channel id: " + str(reference_id)) # Reduce the size of the window for which the correlation will be computed to make space for # the search space around it. channel_window = image[max_shift:-max_shift, max_shift:-max_shift, channel_id] channel_reference = image[:, :, reference_id] # Apply Gaussian blur if requested. if blur_strength: channel_blurred = GaussianBlur(channel_window, (blur_strength, blur_strength), 0) channel_reference_blurred = GaussianBlur(channel_reference, (blur_strength, blur_strength), 0) else: channel_blurred = channel_window channel_reference_blurred = channel_reference # Compute the normalized cross correlation. result = matchTemplate((channel_reference_blurred).astype(float32), channel_blurred.astype(float32), TM_CCORR_NORMED) # Determine the position of the maximum correlation, and compute the corresponding shifts. minVal, maxVal, minLoc, maxLoc = minMaxLoc(result) shift_y = max_shift - maxLoc[1] shift_x = max_shift - maxLoc[0] if abs(shift_y) != max_shift and abs(shift_x) != max_shift: return (shift_y, shift_x) else: raise Error ("No matching shift found within given range") @staticmethod def insert_cross(frame, y_center, x_center, cross_half_len, color): """ Insert a colored cross into an image at a given location. :param frame: Image into which the cross is to be inserted :param y_center: y pixel coordinate of the center of the cross :param x_center: x pixel coordinate of the center of the cross :param cross_half_len: Extension of the cross from center in all four directions :param color: Color, one out of "white", "red", "green" and "blue" :return: - """ shape_y, shape_x = frame.shape[0:2] if color == 'white': rgb = [255, 255, 255] elif color == 'red': rgb = [255, 0, 0] elif color == 'green': rgb = [0, 255, 0] elif color == 'blue': rgb = [0, 0, 255] elif color == 'cyan': rgb = [0, 255, 255] else: rgb = [255, 255, 255] # Be sure not to draw beyond frame boundaries. if 0 <= x_center < shape_x: for y in range(y_center - cross_half_len, y_center + cross_half_len + 1): if 0 <= y < shape_y: frame[y, x_center] = rgb if 0 <= y_center < shape_y: for x in range(x_center - cross_half_len, x_center + cross_half_len + 1): if 0 <= x < shape_x: frame[y_center, x] = rgb @staticmethod def compose_image(image_list, scale_factor=1, border=5): """ Arrange a list of monochrome images horizontally in a single image, with constant gaps in between. If images are of different size, center-align them vertically. Optionally, the resulting image can be re-scaled with the same factor in x and y directions. :param image_list: List containing images (all of the same dtype) :param scale_factor: Scaling factor for the resulting composite image :param border: Width of black border and gaps between images :return: Composite image of the same dtype as the image_list items """ shapes_y = [image.shape[0] for image in image_list] shapes_x = [image.shape[1] for image in image_list] max_shape_y = max(shapes_y) sum_shape_x = sum(shapes_x) # Check if all images of the list are of the same type. type = image_list[0].dtype for image in image_list[1:]: if image.dtype != type: raise ArgumentError("Trying to compose images of different types") # Compute the dimensions of the composite image. The border and the gaps between images are # 5 pixels wide. composite_dim_y = max_shape_y + 2 * border composite_dim_x = sum_shape_x + border * (len(image_list) + 1) # Allocate the composite image. composite = full((composite_dim_y, composite_dim_x), 0, dtype=type) # Copy the images into the composite image. x_pos = border for index, image in enumerate(image_list): y_pos = int((composite_dim_y - shapes_y[index]) / 2) composite[y_pos: y_pos + shapes_y[index], x_pos: x_pos + shapes_x[index]] = image x_pos += shapes_x[index] + border # Return the resized composite image. return resize(composite, None, fx=float(scale_factor), fy=float(scale_factor)) @staticmethod def circle_around(y, x, r): """ Create an iterator which returns y, x pixel coordinates around a given position in an image. Successive elements describe a circle around y, x with distance r (in the 1-norm). The iterator ends when the full circle is traversed. :param x: x coordinate of the center location :param y: y coordinate of the center location :param r: distance of the iterator items from position (y, x) (1-norm, in pixels) :return: - """ # Special case 0: The only iterator element is the center point itself. if r == 0: yield (y, x) # For the general case, set the first element and circle around (y, x) counter-clockwise. j, i = y - r, x - r while i < x + r: i += 1 yield (j, i) while j < y + r: j += 1 yield (j, i) while i > x - r: i -= 1 yield (j, i) while j > y - r: j -= 1 yield (j, i) @staticmethod def write_video(name, frame_list, annotations, fps, depth=8): """ Create a video file from a list of frames. :param name: File name of the video output :param frame_list: List of frames to be written :param annotations: List of text strings to be written into the corresponding frame :param fps: Frames per second of video :param depth: Bit depth of the image "frame", either 8 or 16. :return: - """ # Delete the output file if it exists. try: unlink(name) except: pass # Compute video frame size. frame_height = frame_list[0].shape[0] frame_width = frame_list[0].shape[1] # Define the codec and create VideoWriter object fourcc = VideoWriter_fourcc('D', 'I', 'V', 'X') out = VideoWriter(name, fourcc, fps, (frame_width, frame_height)) # Define font for annotations. font = FONT_HERSHEY_SIMPLEX bottomLeftCornerOfText = (20, 50) fontScale = 1 fontColor = (255, 255, 255) fontThickness = 1 lineType = LINE_AA # For each frame: If monochrome, convert it to three-channel color mode. insert annotation # and write the frame. for index, frame in enumerate(frame_list): # If the frames are 16bit, convert them to 8bit. if depth == 8: frame_8b = frame else: frame_8b = (frame / 255.).astype(uint8) if len(frame.shape) == 2: rgb_frame = stack((frame_8b,) * 3, -1) else: rgb_frame = frame_8b putText(rgb_frame, annotations[index], bottomLeftCornerOfText, font, fontScale, fontColor, fontThickness, lineType) out.write(rgb_frame) out.release() @staticmethod def post_process(image, layers): """ Apply all postprocessing layers to the input image "image". If the image is in color mode, the postprocessing is computed either in BGR mode or on the luminance channel only. All computations are performed in 32bit floating point mode. :param image: Input image, either BGR or Grayscale (16bit uint). :param layers: postprocessing layers with all parameters. :return: Processed image in the same 16bit uint format as the input image. """ # Check if the original image is selected (version 0). In this case nothing is to be done. if not layers: return image.astype(uint16) # Convert the image to 32bit floating point format. input_image = image.astype(float32) color = len(image.shape) == 3 # If the luminance channel is to be processed only, extract the luminance channel. if color and layers[0].luminance_only: input_image_hsv = cvtColor(input_image, COLOR_BGR2HSV) layer_input = input_image_hsv[:, :, 2] else: layer_input = input_image # Go through all layers and apply the sharpening filters. for layer_index, layer in enumerate(layers): # Bilateral filter is needed: if abs(layer.bi_fraction) > 1.e-5: layer_bilateral = bilateralFilter(layer_input, 0, layer.bi_range * 256., layer.radius / 3., borderType=BORDER_DEFAULT) # Gaussian filter is needed: if abs(layer.bi_fraction - 1.) > 1.e-5: layer_gauss = GaussianBlur(layer_input, (0, 0), layer.radius / 3., borderType=BORDER_DEFAULT) # Compute the input for the next layer. First case: bilateral only. if abs(layer.bi_fraction - 1.) <= 1.e-5: next_layer_input = layer_bilateral # Case Gaussian only. elif abs(layer.bi_fraction) <= 1.e-5: next_layer_input = layer_gauss # Mixed case. else: next_layer_input = layer_bilateral * layer.bi_fraction + \ layer_gauss * (1. - layer.bi_fraction) layer_component_before_denoise = layer_input - next_layer_input # If denoising is chosen for this layer, apply a Gaussian filter. if layer.denoise > 1.e-5: layer_component = GaussianBlur(layer_component_before_denoise, (0, 0), layer.radius / 3., borderType=BORDER_DEFAULT) * layer.denoise + \ layer_component_before_denoise * (1. - layer.denoise) else: layer_component = layer_component_before_denoise # Accumulate the contributions from all layers. On the first layer initialize the # summation buffer. if layer_index: components_accumulated += layer_component * layer.amount else: components_accumulated = layer_component * layer.amount layer_input = next_layer_input # After all layers are accumulated, finally add the maximally blurred image. components_accumulated += next_layer_input # Reduce the value range so that they fit into 16bit uint, and convert to uint16. In case # the luminance channel was processed only, insert the processed luminance channel into the # HSV representation of the original image, and convert back to BGR. if color and layers[0].luminance_only: input_image_hsv[:, :, 2] = components_accumulated return cvtColor(input_image_hsv.clip(min=0., max=65535.), COLOR_HSV2BGR).astype(uint16) else: return components_accumulated.clip(min=0., max=65535.).astype(uint16) @staticmethod def gaussian_sharpen(input_image, amount, radius, luminance_only=False): """ Sharpen an image with a Gaussian kernel. The input image can be B/W or color. :param input_image: Input image, type uint16 :param amount: Amount of sharpening :param radius: Radius of Gaussian kernel (in pixels) :param luminance_only: True, if only the luminance channel of a color image is to be sharpened. Default is False. :return: The sharpened image (B/W or color, as input), type uint16 """ color = len(input_image.shape) == 3 # Translate the kernel radius into standard deviation. sigma = radius / 3 # Convert the image to floating point format. image = input_image.astype(float32) # Special case: Only sharpen the luminance channel of a color image. if color and luminance_only: hsv = cvtColor(image, COLOR_BGR2HSV) luminance = hsv[:, :, 2] # Apply a Gaussian blur filter, subtract it from the original image, and add a multiple # of this correction to the original image. Clip values out of range. luminance_blurred = GaussianBlur(luminance, (0, 0), sigma, borderType=BORDER_DEFAULT) hsv[:, :, 2] = (luminance + amount * (luminance - luminance_blurred)).clip(min=0., max=65535.) # Convert the image back to uint16. return cvtColor(hsv, COLOR_HSV2BGR).astype(uint16) # General case: Treat the entire image (B/W or color 16bit mode). else: image_blurred = GaussianBlur(image, (0, 0), sigma, borderType=BORDER_DEFAULT) return (image + amount * (image - image_blurred)).clip(min=0., max=65535.).astype( uint16) @staticmethod def gaussian_blur(input_image, amount, radius, luminance_only=False): """ Soften an image with a Gaussian kernel. The input image can be B/W or color. :param input_image: Input image, type uint16 :param amount: Amount of blurring, between 0. and 1. :param radius: Radius of Gaussian kernel (in pixels) :param luminance_only: True, if only the luminance channel of a color image is to be blurred. Default is False. :return: The blurred image (B/W or color, as input), type uint16 """ color = len(input_image.shape) == 3 # Translate the kernel radius into standard deviation. sigma = radius / 3 # Convert the image to floating point format. image = input_image.astype(float32) # Special case: Only blur the luminance channel of a color image. if color and luminance_only: hsv = cvtColor(image, COLOR_BGR2HSV) luminance = hsv[:, :, 2] # Apply a Gaussian blur filter, subtract it from the original image, and add a multiple # of this correction to the original image. Clip values out of range. luminance_blurred = GaussianBlur(luminance, (0, 0), sigma, borderType=BORDER_DEFAULT) hsv[:, :, 2] = (luminance_blurred*amount + luminance*(1.-amount)).clip(min=0., max=65535.) # Convert the image back to uint16. return cvtColor(hsv, COLOR_HSV2BGR).astype(uint16) # General case: Treat the entire image (B/W or color 16bit mode). else: image_blurred = GaussianBlur(image, (0, 0), sigma, borderType=BORDER_DEFAULT) return (image_blurred*amount + image*(1.-amount)).clip(min=0., max=65535.).astype( uint16) @staticmethod def wavelet_sharpen(input_image, amount, radius): """ Sharpen a B/W or color image with wavelets. The underlying algorithm was taken from the Gimp wavelet plugin, originally written in C and published under the GPLv2+ license at: https://github.com/mrossini-ethz/gimp-wavelet-sharpen/blob/master/src/wavelet.c :param input_image: Input image (B/W or color), type uint16 :param amount: Amount of sharpening :param radius: Radius in pixels :return: Sharpened image, same format as input image """ height, width = input_image.shape[:2] color = len(input_image.shape) == 3 # Allocate workspace: Three complete images, plus 1D object with length max(row, column). if color: fimg = empty((3, height, width, 3), dtype=float32) temp = zeros((max(width, height), 3), dtype=float32) else: fimg = empty((3, height, width), dtype=float32) temp = zeros(max(width, height), dtype=float32) # Convert input image to floats. fimg[0] = input_image / 65535 # Start with level 0. Store its Laplacian on level 1. The operator is separated in a # column and a row operator. hpass = 0 for lev in range(5): # Highpass and lowpass levels use image indices 1 and 2 in alternating mode to save # space. lpass = ((lev & 1) + 1) if color: for row in range(height): Miscellaneous.mexican_hat_color(temp, fimg[hpass][row, :, :], width, 1 << lev) fimg[lpass][row, :, :] = temp[:width, :] * 0.25 for col in range(width): Miscellaneous.mexican_hat_color(temp, fimg[lpass][:, col, :], height, 1 << lev) fimg[lpass][:, col, :] = temp[:height, :] * 0.25 else: for row in range(height): Miscellaneous.mexican_hat(temp, fimg[hpass][row, :], width, 1 << lev) fimg[lpass][row, :] = temp[:width] * 0.25 for col in range(width): Miscellaneous.mexican_hat(temp, fimg[lpass][:, col], height, 1 << lev) fimg[lpass][:, col] = temp[:height] * 0.25 # Compute the amount of the correction at the current level. amt = amount * exp(-(lev - radius) * (lev - radius) / 1.5) + 1. fimg[hpass] -= fimg[lpass] fimg[hpass] *= amt # Accumulate all corrections in the first workspace image. if hpass: fimg[0] += fimg[hpass] hpass = lpass # At the end add the coarsest level and convert back to 16bit integer format. fimg[0] = ((fimg[0] + fimg[lpass]) * 65535.).clip(min=0., max=65535.) return fimg[0].astype(uint16) @staticmethod def mexican_hat(temp, base, size, sc): """ Apply a 1D strided second derivative to a row or column of a B/W image. Store the result in the temporary workspace "temp". :param temp: Workspace (type float32), length at least "size" elements :param base: Input image (B/W), Type float32 :param size: Length of image row / column :param sc: Stride (power of 2) of operator :return: - """ # Special case at begin of row/column. Full operator not applicable. temp[:sc] = 2 * base[:sc] + base[sc:0:-1] + base[sc:2 * sc] # Apply the full operator. temp[sc:size - sc] = 2 * base[sc:size - sc] + base[:size - 2 * sc] + base[2 * sc:size] # Special case at end of row/column. The full operator is not applicable. temp[size - sc:size] = 2 * base[size - sc:size] + base[size - 2 * sc:size - sc] + \ base[size - 2:size - 2 - sc:-1] @staticmethod def mexican_hat_color(temp, base, size, sc): """ Apply a 1D strided second derivative to a row or column of a color image. Store the result in the temporary workspace "temp". :param temp: Workspace (type float32), length at least "size" elements (first dimension) times 3 colors (second dimension). :param base: Input image (color), Type float32 :param size: Length of image row / column :param sc: Stride (power of 2) of operator :return: - """ # Special case at begin of row/column. Full operator not applicable. temp[:sc, :] = 2 * base[:sc, :] + base[sc:0:-1, :] + base[sc:2 * sc, :] # Apply the full operator. temp[sc:size - sc, :] = 2 * base[sc:size - sc, :] + base[:size - 2 * sc, :] + base[ 2 * sc:size, :] # Special case at end of row/column. The full operator is not applicable. temp[size - sc:size, :] = 2 * base[size - sc:size, :] + base[size - 2 * sc:size - sc, :] + \ base[size - 2:size - 2 - sc:-1, :] @staticmethod def protocol(string, logfile, precede_with_timestamp=True): """ Print a message (optionally plus time stamp) to standard output. If it is requested to store a logfile with the stacked image, write the string to that file as well. :param string: Message to be printed after the time stamp :param logfile: logfile or None (no logging) :return: - """ # Precede the text with a time stamp and print it to stdout. Note that stdout may be # redirected to a file. if precede_with_timestamp: output_string = '{0} {1}'.format(datetime.now().strftime("%H-%M-%S.%f")[:-5], string) else: output_string = string print(output_string) stdout.flush() # If a logfile per stacked image was requested, write the string to that file as well. if logfile: logfile.write(output_string + "\n") @staticmethod def print_stacking_parameters(configuration, logfile): """ Print a table with all stacking parameters. :param configuration: Object holding all configuration parameters. :param logfile: logfile or None (no logging) :return: - """ # Compile all parameters and their values to be printed. if configuration.global_parameters_maximum_memory_active: parameters = [["RAM limit set explicitly (GBytes)", str(configuration.global_parameters_maximum_memory_amount)]] else: parameters = [["RAM limit set explicitly", "False"]] if configuration.global_parameters_buffering_level == -1: parameters += [["Buffering level", "auto"]] else: parameters += [["Buffering level set explicitly", str(configuration.global_parameters_buffering_level)]] parameters += [["Debayering default", configuration.frames_debayering_default], ["Debayering method", configuration.frames_debayering_method], ["Noise level (add Gaussian blur)", str(configuration.frames_gauss_width)], ["Frame stabilization mode", configuration.align_frames_mode]] # The following parameters are only active in "Surface" mode. if configuration.align_frames_mode == "Surface": parameters += [ ["Automatic frame stabilization", str(configuration.align_frames_automation)], ["Stabilization patch size (% of frame size)", str(int(round(100. / configuration.align_frames_rectangle_scale_factor)))], ["Stabilization search width (pixels)", str(configuration.align_frames_search_width)]] # Continue with general parameters. parameters += [["Percentage of best frames for reference frame computation", str(configuration.align_frames_average_frame_percent)], ["Object is changing fast (e.g. Jupiter, Sun)", str(configuration.align_frames_fast_changing_object)], ["Alignment box width (pixels)", str(2 * configuration.alignment_points_half_box_width)], ["Max. alignment search width (pixels)", str(configuration.alignment_points_search_width)], ["Minimum structure", str(configuration.alignment_points_structure_threshold)], ["Minimum brightness", str(configuration.alignment_points_brightness_threshold)]] # The stack size is defined either as frame number or percentage. if configuration.alignment_points_frame_percent > 0: parameters += [["Percentage of best frames to be stacked", str(configuration.alignment_points_frame_percent)], ["Normalize frame brightness", str(configuration.frames_normalization)]] else: parameters += [["Number of best frames to be stacked", str(configuration.alignment_points_frame_number)], ["Normalize frame brightness", str(configuration.frames_normalization)]] # If brightness normalization is checked, add the black cut-off value. if configuration.frames_normalization: parameters += [ ["Normalization black cut-off", str(configuration.frames_normalization_threshold)]] # If drizzling is active, add the factor. if configuration.stack_frames_drizzle_factor_string != "Off": parameters += [["Drizzle factor in stacking", str(configuration.stack_frames_drizzle_factor_string)]] output_string = "\n Stacking parameters: " \ " | Value |\n" \ " " \ "---------------------------------------------------------------------------------------------" \ "\n " # Extend the output string with a line for every parameter to be printed. for line in parameters: output_string += " {0:60s} | {1:29s}|\n ".format(line[0], line[1]) # Write the complete table. Miscellaneous.protocol(output_string, logfile, precede_with_timestamp=False) @staticmethod def print_postproc_parameters(postproc_version, logfile): """ Print a table with postprocessing layer info for the selected postprocessing version. :param postproc_version: Object holding postprocessing parameters for the selected version. :param logfile: logfile or None (no logging) :return: - """ # Test if an RGB correction has been applied. if postproc_version.shift_red != (0., 0.) or postproc_version.shift_blue != (0., 0.): # Find out if the RGB correction was done automatically. if postproc_version.rgb_automatic: intro = " Automatic RGB correction, " else: intro = " Manual RGB correction, " (shift_red_y, shift_red_x) = postproc_version.shift_red (shift_blue_y, shift_blue_x) = postproc_version.shift_blue n_digits = [0, 1, 2][postproc_version.rgb_resolution_index] if shift_red_y >= 0.: dir_red_y = " pixels down" else: dir_red_y = " pixels up" if shift_red_x >= 0.: dir_red_x = " pixels right" else: dir_red_x = " pixels left" if shift_blue_y >= 0.: dir_blue_y = " pixels down" else: dir_blue_y = " pixels up" if shift_blue_x >= 0.: dir_blue_x = " pixels right" else: dir_blue_x = " pixels left" # Special case 0 digits: In this case the number of digits must be omitted. # If the round function is called with "n_digits=0", the result still has one digit # after the decimal point. if n_digits: Miscellaneous.protocol( intro + "red channel shifted " + str(round(abs(shift_red_y), n_digits)) + dir_red_y + ", " + str(round(abs(shift_red_x), n_digits)) + dir_red_x + ", blue channel shifted " + str(round(abs(shift_blue_y), n_digits)) + dir_blue_y + ", " + str(round(abs(shift_blue_x), n_digits)) + dir_blue_x + ".", logfile, precede_with_timestamp=False) else: Miscellaneous.protocol( intro + "red channel shifted " + str(round(abs(shift_red_y))) + dir_red_y + ", " + str(round(abs(shift_red_x))) + dir_red_x + ", blue channel shifted " + str(round(abs(shift_blue_y))) + dir_blue_y + ", " + str(round(abs(shift_blue_x))) + dir_blue_x + ".", logfile, precede_with_timestamp=False) if postproc_version.layers: output_string = " Postprocessing method: " + postproc_version.layers[0].postproc_method + "\n\n" + \ " Layer | Radius | Amount | Bi fraction (%) | Bi range | Denoise (%) | Luminance only |\n" \ " ------------------------------------------------------------------------------------------------------------------" \ "\n " # Extend the three table lines up to the max index. for index, layer in enumerate(postproc_version.layers): output_string += " {0:3d} | {1:5.2f} | {2:6.2f} | {3:4.0f} | {4:5.1f} | {5:4.0f} | {6:8s} |" \ "\n ".format(index + 1, layer.radius, layer.amount, layer.bi_fraction*100., layer.bi_range, layer.denoise*100., str(layer.luminance_only)) Miscellaneous.protocol(output_string, logfile, precede_with_timestamp=False) if __name__ == "__main__": # Test program for routine search_local_match # Set the size of the image frame. frame_height = 960 frame_width = 1280 # Initialize the frame with a wave-like pattern in x and y directions. x_max = 30. x_vec = arange(0., x_max, x_max / frame_width) y_max = 25. y_vec = arange(0., y_max, y_max / frame_height) frame = empty((frame_height, frame_width)) for y_j, y in enumerate(y_vec): for x_i, x in enumerate(x_vec): frame[y_j, x_i] = sin(y) * sin(x) # Set the size and location of the reference frame window and cut it out from the frame. window_height = 40 window_width = 50 reference_y_low = 400 reference_x_low = 500 reference_box = frame[reference_y_low:reference_y_low + window_height, reference_x_low:reference_x_low + window_width] # Set the true displacement vector to be checked against the result of the search function. displacement_y = 1 displacement_x = -2 # The start point for the local search is offset from the true matching point. y_low = reference_y_low + displacement_y y_high = y_low + window_height x_low = reference_x_low + displacement_x x_high = x_low + window_width # Set the radius of the search area. search_width = 20 sampling_stride = 1 dev_table = empty((2 * search_width + 1, 2 * search_width + 1), dtype=float32) # First try the standard method "search_local_match". Compute the displacement vector, and # print a comparison of the true and computed values. start = time() rep_count = 1 for iter in range(rep_count): [dy, dx], dev_r = Miscellaneous.search_local_match(reference_box, frame, y_low, y_high, x_low, x_high, search_width, sampling_stride, sub_pixel=False) end = time() print("Standard method, true displacements: " + str([displacement_y, displacement_x]) + ", computed: " + str([dy, dx]) + ", execution time (s): " + str((end - start) / rep_count)) # As a comparison, do the same using the steepest descent method. start = time() for iter in range(rep_count): [dy, dx], dev_r = Miscellaneous.search_local_match_gradient(reference_box, frame, y_low, y_high, x_low, x_high, search_width, sampling_stride, dev_table) end = time() print("steepest descent, true displacements: " + str([displacement_y, displacement_x]) + ", computed: " + str([dy, dx]) + ", execution time (s): " + str((end - start) / rep_count)) # Now test the alternative method in comparison to the straight-forward one. First, # compute the reference frame box stack. The initialization has to be performed only once for # each alignment point. It can be reused for all frames containing the AP. start = time() for iter in range(rep_count): reference_stack, displacements, radius_start = \ Miscellaneous.search_local_match_init(frame, reference_y_low, reference_y_low + window_height, reference_x_low, reference_x_low + window_width, search_width) end = time() print("Match initialization time (s): " + str((end - start) / rep_count)) # Compute the displacement. The result should be the same as above. start = time() for iter in range(rep_count): [dy, dx], dev_r = Miscellaneous.search_local_match_execute( frame[y_low:y_high, x_low:x_high], reference_stack, displacements, radius_start) end = time() print("Match execution, true displacements: " + str( [displacement_y, displacement_x]) + ", computed: " + str( [dy, dx]) + ", execution time (s): " + str((end - start) / rep_count))