satdump/src-core/image/processing.cpp

610 lines
No EOL
26 KiB
C++

#include "image_utils.h"
#include <algorithm>
#include <cfloat>
#include <cmath>
#include <omp.h>
#include <thread>
#include "logger.h"
#include "utils/time.h"
#include "core/opencl.h"
#include "core/resources.h"
namespace satdump
{
namespace image
{
int percentile(int *array, int size, float percentile)
{
float number_percent = (size + 1) * percentile / 100.0f;
if (number_percent <= 1) // ==
return array[0];
else if (number_percent >= size) // ==
return array[size - 1];
else
return array[(int)number_percent - 1] + (number_percent - (int)number_percent) * (array[(int)number_percent] - array[(int)number_percent - 1]);
}
inline int wraparound(int metric, int attempt)
{
if (attempt < 0)
attempt += metric;
if (metric <= attempt)
attempt -= metric;
return attempt;
}
void white_balance(Image &img, float percentileValue)
{
const float maxVal = img.maxval();
const size_t d_height = img.height();
const size_t d_width = img.width();
int *sorted_array = new int[d_height * d_width];
for (int c = 0; c < img.channels(); c++)
{
// Load the whole image band into our array
for (size_t i = 0; i < d_height * d_width; i++)
sorted_array[i] = img.get(c, i);
// Sort it
std::sort(&sorted_array[0], &sorted_array[d_width * d_height]);
// Get percentiles
int percentile1 = percentile(sorted_array, d_width * d_height, percentileValue);
int percentile2 = percentile(sorted_array, d_width * d_height, 100.0f - percentileValue);
for (size_t i = 0; i < d_width * d_height; i++)
{
long balanced = (img.get(c, i) - percentile1) * maxVal / (percentile2 - percentile1);
if (balanced < 0)
balanced = 0;
else if (balanced > maxVal)
balanced = maxVal;
img.set(c, i, balanced);
}
}
delete[] sorted_array;
}
void median_blur(Image &img)
{
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
Image orig = img;
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
int values[5];
size_t idx = (size_t)x * w + y;
values[0] = values[1] = values[2] = values[3] = values[4] = orig.get(c, idx);
if (x != 0) values[1] = orig.get(c, (size_t)(x - 1) * w + y);
if (y != 0) values[2] = orig.get(c, (size_t)x * w + (y - 1));
if (x != (int64_t)h - 1) values[3] = orig.get(c, (size_t)(x + 1) * w + y);
if (y != w - 1) values[4] = orig.get(c, (size_t)x * w + (y + 1));
std::nth_element(values, values + 2, values + 5);
img.set(c, idx, values[2]);
}
}
}
}
void kuwahara_filter(Image &img)
{
const int radius = 1;
const float num_pixels = (float)((radius + 1) * (radius + 1));
const int d_channels = img.channels();
const size_t d_width = img.width();
const size_t d_height = img.height();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
Image orig = img;
for (int c = 0; c < d_channels; c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t y = 0; y < (int64_t)d_height; y++)
{
for (size_t x = 0; x < d_width; x++)
{
float average[4] = {0};
float variance[4] = {0};
// Calculate values for the four regions
for (int j = -radius; j <= 0; ++j)
for (int i = -radius; i <= 0; ++i)
average[0] += orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j));
average[0] /= num_pixels;
for (int j = -radius; j <= 0; ++j)
for (int i = -radius; i <= 0; ++i)
variance[0] += pow(orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j)) - average[0], 2);
for (int j = -radius; j <= 0; ++j)
for (int i = 0; i <= radius; ++i)
average[1] += orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j));
average[1] /= num_pixels;
for (int j = -radius; j <= 0; ++j)
for (int i = 0; i <= radius; ++i)
variance[1] += pow(orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j)) - average[1], 2);
for (int j = 0; j <= radius; ++j)
for (int i = 0; i <= radius; ++i)
average[2] += orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j));
average[2] /= num_pixels;
for (int j = 0; j <= radius; ++j)
for (int i = 0; i <= radius; ++i)
variance[2] += pow(orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j)) - average[2], 2);
for (int j = 0; j <= radius; ++j)
for (int i = -radius; i <= 0; ++i)
average[3] += orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j));
average[3] /= num_pixels;
for (int j = 0; j <= radius; ++j)
for (int i = -radius; i <= 0; ++i)
variance[3] += pow(orig.get(c, wraparound(d_width, x + i), wraparound(d_height, y + j)) - average[3], 2);
// Find the region with the smallest variance and use its mean as the new pixel value
float min_sigma2 = FLT_MAX;
for (int k = 0; k < 4; k++)
{
variance[k] /= num_pixels - 1;
if (variance[k] < 0) variance[k] = -variance[k];
if (variance[k] < min_sigma2)
{
min_sigma2 = variance[k];
img.set(c, x, y, (int)(average[k] + 0.5f));
}
}
}
}
}
}
void equalize(Image &img, bool per_channel)
{
for (int c = 0; c < (per_channel ? img.channels() : 1); c++)
{
if (c == 3) // Do not individual equalize alpha channel
break;
const int nlevels = img.maxval() + 1;
size_t size = img.width() * img.height() * (per_channel ? 1 : img.channels());
// Init histogram buffer
int *histogram = new int[nlevels];
for (int i = 0; i < nlevels; i++)
histogram[i] = 0;
// Compute histogram
for (size_t px = 0; px < size; px++)
histogram[img.get(c, px)]++;
// Cummulative histogram
int *cummulative_histogram = new int[nlevels];
cummulative_histogram[0] = histogram[0];
for (int i = 1; i < nlevels; i++)
cummulative_histogram[i] = histogram[i] + cummulative_histogram[i - 1];
// Scaling
int *scaling = new int[nlevels];
for (int i = 0; i < nlevels; i++)
scaling[i] = round(cummulative_histogram[i] * (float(nlevels - 1) / size));
// Apply
for (size_t px = 0; px < size; px++)
img.set(c, px, img.clamp(scaling[img.get(c, px)]));
// Cleanup
delete[] cummulative_histogram;
delete[] scaling;
delete[] histogram;
}
}
void normalize(Image &img)
{
int max = 0;
int min = img.maxval();
// Get min and max
for (size_t i = 0; i < img.size(); i++)
{
int val = img.get(i);
if (val > max)
max = val;
if (val < min)
min = val;
}
if (abs(max - min) == 0) // Avoid division by 0
return;
// Compute scaling factor
float factor = img.maxval() / float(max - min);
// Scale entire image
for (size_t i = 0; i < img.size(); i++)
img.set(i, img.clamp((img.get(i) - min) * factor));
}
void linear_invert(Image &img)
{
for (size_t i = 0; i < img.size(); i++)
img.set(i, img.maxval() - img.get(i));
}
void simple_despeckle(Image &img, int thresold)
{
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
Image orig = img;
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
unsigned short current = orig.get(c, (size_t)x * w + y);
unsigned short below = x + 1 == (int64_t)h ? 0 : orig.get(c, (size_t)(x + 1) * w + y);
unsigned short left = y - 1 == -1 ? 0 : orig.get(c, (size_t)x * w + (y - 1));
unsigned short right = y + 1 == w ? 0 : orig.get(c, (size_t)x * w + (y + 1));
if ((current - left > thresold && current - right > thresold) || (current - below > thresold && current - right > thresold))
{
img.set(c, (size_t)x * w + y, (right + left) / 2);
}
}
}
}
}
void switching_median(Image &img, int threshold)
{
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
int values[5];
size_t idx = (size_t)x * w + y;
values[0] = values[1] = values[2] = values[3] = values[4] = orig.get(c, idx);
if (x != 0) values[1] = orig.get(c, (size_t)(x - 1) * w + y);
if (y != 0) values[2] = orig.get(c, (size_t)x * w + (y - 1));
if (x != (int64_t)h - 1) values[3] = orig.get(c, (size_t)(x + 1) * w + y);
if (y != w - 1) values[4] = orig.get(c, (size_t)x * w + (y + 1));
std::nth_element(values, values + 2, values + 5);
int current = orig.get(c, idx);
int median = values[2];
if (abs(current - median) > threshold)
img.set(c, idx, median);
}
}
}
}
void adaptive_median(Image &img, int threshold_strong)
{
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
int v3[9];
int count3 = 0;
for (int i = -1; i <= 1; i++)
for (int j = -1; j <= 1; j++)
if (x + i >= 0 && x + i < (int64_t)h && y + j >= 0 && y + j < w)
v3[count3++] = orig.get(c, (size_t)(x + i) * w + (y + j));
std::nth_element(v3, v3 + count3 / 2, v3 + count3);
int m3 = v3[count3 / 2];
int current = orig.get(c, (size_t)x * w + y);
if (abs(current - m3) > threshold_strong)
{
int v5[25];
int count5 = 0;
for (int i = -2; i <= 2; i++)
for (int j = -2; j <= 2; j++)
if (x + i >= 0 && x + i < (int64_t)h && y + j >= 0 && y + j < w)
v5[count5++] = orig.get(c, (size_t)(x + i) * w + (y + j));
std::nth_element(v5, v5 + count5 / 2, v5 + count5);
img.set(c, (size_t)x * w + y, v5[count5 / 2]);
}
else
{
img.set(c, (size_t)x * w + y, m3);
}
}
}
}
}
void bilateral_filter(Image &img, int radius, float sigma_intensity)
{
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
float sigma_space = radius / 2.0f;
// Automatically scale sigma_intensity to the image dynamic range
// (Assuming sigma_intensity 25 is "medium" for 8-bit images)
float scaled_sigma = sigma_intensity * (img.maxval() / 255.0f);
if (scaled_sigma < 0.1f) scaled_sigma = 0.1f;
// Precalculate spatial weights
std::vector<float> spatial_weights((2 * radius + 1) * (2 * radius + 1));
for (int i = -radius; i <= radius; i++)
for (int j = -radius; j <= radius; j++)
spatial_weights[(i + radius) * (2 * radius + 1) + (j + radius)] = exp(-(float)(i * i + j * j) / (2 * sigma_space * sigma_space));
// Intensity difference LUT (Up to 16-bit range)
std::vector<float> range_lut(65536);
double inv_sigma_range = 1.0 / (2.0 * (double)scaled_sigma * (double)scaled_sigma);
for (int i = 0; i < 65536; i++)
range_lut[i] = exp(-((double)i * (double)i) * inv_sigma_range);
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
float sum_weights = 0;
float sum_values = 0;
int current = orig.get(c, (size_t)x * w + y);
for (int i = -radius; i <= radius; i++)
{
for (int j = -radius; j <= radius; j++)
{
int64_t nx = x + i;
int64_t ny = (int64_t)y + j;
if (nx >= 0 && nx < (int64_t)h && ny >= 0 && ny < (int64_t)w)
{
int neighbor = orig.get(c, (size_t)nx * w + (size_t)ny);
int diff = abs(current - neighbor);
float w_space = spatial_weights[(i + radius) * (2 * radius + 1) + (j + radius)];
float w_range = (diff < 65536) ? range_lut[diff] : 0.0f;
float weight = w_space * w_range;
sum_weights += weight;
sum_values += neighbor * weight;
}
}
}
if (sum_weights > 0.00001f)
img.set(c, (size_t)x * w + y, (int)(sum_values / sum_weights + 0.5f));
}
}
}
}
void simple_inpainting(Image &img, int threshold)
{
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
for (size_t y = 0; y < w; y++)
{
int v3[9];
int count3 = 0;
for (int i = -1; i <= 1; i++)
for (int j = -1; j <= 1; j++)
if (x + i >= 0 && x + i < (int64_t)h && y + j >= 0 && y + j < w)
v3[count3++] = orig.get(c, (size_t)(x + i) * w + (y + j));
std::nth_element(v3, v3 + count3 / 2, v3 + count3);
int m3 = v3[count3 / 2];
int current = orig.get(c, (size_t)x * w + y);
if (abs(current - m3) > threshold)
{
float sum = 0;
int count = 0;
for (int i = -1; i <= 1; i++)
{
for (int j = -1; j <= 1; j++)
{
if (i == 0 && j == 0) continue;
int64_t nx = x + i;
int64_t ny = (int64_t)y + j;
if (nx >= 0 && nx < (int64_t)h && ny >= 0 && ny < (int64_t)w)
{
int val = orig.get(c, (size_t)nx * w + (size_t)ny);
if (abs(val - m3) <= threshold)
{
sum += val;
count++;
}
}
}
}
if (count > 0)
img.set(c, (size_t)x * w + y, (int)(sum / count + 0.5f));
else
img.set(c, (size_t)x * w + y, m3);
}
else
{
img.set(c, (size_t)x * w + y, current);
}
}
}
}
}
void selective_impulse_filter(Image &img, int threshold, int window_size)
{
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
int r = (window_size - 1) / 2;
for (int c = 0; c < img.channels(); c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t x = 0; x < (int64_t)h; x++)
{
std::vector<int> values;
values.reserve(window_size * window_size);
for (size_t y = 0; y < w; y++)
{
values.clear();
int current = orig.get(c, (size_t)x * w + y);
for (int i = -r; i <= r; i++)
{
for (int j = -r; j <= r; j++)
{
int64_t nx = x + i;
int64_t ny = (int64_t)y + j;
if (nx >= 0 && nx < (int64_t)h && ny >= 0 && ny < (int64_t)w)
{
values.push_back(orig.get(c, (size_t)nx * w + (size_t)ny));
}
}
}
// We need the entire sorted window to find gaps
std::sort(values.begin(), values.end());
int median = values[values.size() / 2];
// A pixel is noise ONLY if it's strictly an isolated outlier
bool updated = false;
// Check if current is the absolute minimum AND there's a big gap to the next value
if (current == values[0])
{
// If it's a lone black pixel, the NEXT value in sorted list will be much higher
int next_nearest = values[1];
if ((next_nearest - current) > threshold)
{
img.set(c, (size_t)x * w + y, median);
updated = true;
}
}
// Check if current is the absolute maximum AND there's a big gap to the previous value
if (!updated && current == values[values.size() - 1])
{
// If it's a lone white pixel, the PREVIOUS value in sorted list will be much lower
int prev_nearest = values[values.size() - 2];
if ((current - prev_nearest) > threshold)
{
img.set(c, (size_t)x * w + y, median);
}
}
}
}
}
}
void scanline_noise_remover(Image &img, int threshold, int radius)
{
double t0 = getTime();
Image orig = img;
size_t h = img.height();
size_t w = img.width();
int channels = img.channels();
int num_threads = std::max(1, (int)(std::thread::hardware_concurrency() * 0.75));
double sigma_multiplier = threshold / 10.0;
for (int c = 0; c < channels; c++)
{
#pragma omp parallel for num_threads(num_threads)
for (int64_t y = 0; y < (int64_t)h; y++)
{
double sum = 0;
double sum_sq = 0;
int count = 0;
// Helper to get pixel with border clamping
auto get_pix = [&](int64_t ix) {
if (ix < 0) ix = 0;
if (ix >= (int64_t)w) ix = w - 1;
return (double)orig.get(c, (size_t)ix, (size_t)y);
};
// Initial window setup for the start of the scanline
for (int i = -radius; i <= radius; i++)
{
double val = get_pix(i);
sum += val;
sum_sq += val * val;
count++;
}
for (size_t x = 0; x < w; x++)
{
double mean = sum / count;
double variance = (sum_sq / count) - (mean * mean);
if (variance < 0) variance = 0;
double std_dev = sqrt(variance);
if (std_dev < 1.0) std_dev = 1.0;
int current = (int)get_pix(x);
if (abs(current - mean) > sigma_multiplier * std_dev)
{
int previous = (x > 0) ? (int)get_pix(x - 1) : (int)mean;
img.set(c, x, y, previous);
}
// Sliding window update
double old_val = get_pix((int64_t)x - radius);
double new_val = get_pix((int64_t)x + radius + 1);
sum = sum - old_val + new_val;
sum_sq = sum_sq - (old_val * old_val) + (new_val * new_val);
}
}
}
logger->info("Scanline noise removal took: %7.2f ms", (getTime() - t0) * 1000.0);
}
} // namespace image
} // namespace satdump