openmv/lib/imlib/stats.c
2025-07-31 13:13:45 -07:00

1207 lines
54 KiB
C

/*
* SPDX-License-Identifier: MIT
*
* Copyright (C) 2013-2024 OpenMV, LLC.
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
* THE SOFTWARE.
*
* Statistics functions.
*/
#include "imlib.h"
#ifdef IMLIB_ENABLE_GET_SIMILARITY
typedef struct imlib_similarity_line_op_state {
bool dssim;
int *sumBucketsOfX, *sumBucketsOfY, *sum2BucketsOfX, *sum2BucketsOfY, *sum2Buckets;
float similarity_sum, similarity_sum_2, similarity_min, similarity_max;
int lines_processed, lines;
} imlib_similarity_line_op_state_t;
static void imlib_similarity_line_op(int x, int x_end, int y_row, imlib_draw_row_data_t *data) {
imlib_similarity_line_op_state_t *state = data->callback_arg;
float c1 = 0, c2 = 0;
int x_start = x;
switch (data->dst_img->pixfmt) {
case PIXFORMAT_BINARY: {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(data->dst_img, y_row);
uint32_t *other_row_ptr = (uint32_t *) data->dst_row_override;
for (; x < x_end; x++) {
int pixel = IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x);
int other_pixel = IMAGE_GET_BINARY_PIXEL_FAST(other_row_ptr, x);
int bucket = (x - x_start) / 8;
state->sumBucketsOfX[bucket] += pixel;
state->sumBucketsOfY[bucket] += other_pixel;
state->sum2BucketsOfX[bucket] += pixel * pixel;
state->sum2BucketsOfY[bucket] += other_pixel * other_pixel;
state->sum2Buckets[bucket] += pixel * other_pixel;
}
c1 = (COLOR_BINARY_MAX * 0.01f) * (COLOR_BINARY_MAX * 0.01f);
c2 = (COLOR_BINARY_MAX * 0.03f) * (COLOR_BINARY_MAX * 0.03f);
break;
}
case PIXFORMAT_GRAYSCALE: {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(data->dst_img, y_row);
uint8_t *other_row_ptr = (uint8_t *) data->dst_row_override;
for (; x < x_end; x++) {
int pixel = IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x);
int other_pixel = IMAGE_GET_GRAYSCALE_PIXEL_FAST(other_row_ptr, x);
int bucket = (x - x_start) / 8;
state->sumBucketsOfX[bucket] += pixel;
state->sumBucketsOfY[bucket] += other_pixel;
state->sum2BucketsOfX[bucket] += pixel * pixel;
state->sum2BucketsOfY[bucket] += other_pixel * other_pixel;
state->sum2Buckets[bucket] += pixel * other_pixel;
}
c1 = (COLOR_GRAYSCALE_MAX * 0.01f) * (COLOR_GRAYSCALE_MAX * 0.01f);
c2 = (COLOR_GRAYSCALE_MAX * 0.03f) * (COLOR_GRAYSCALE_MAX * 0.03f);
break;
}
case PIXFORMAT_RGB565: {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(data->dst_img, y_row);
uint16_t *other_row_ptr = (uint16_t *) data->dst_row_override;
for (; x < x_end; x++) {
int pixel = COLOR_RGB565_TO_Y(IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x));
int other_pixel = COLOR_RGB565_TO_Y(IMAGE_GET_RGB565_PIXEL_FAST(other_row_ptr, x));
int bucket = (x - x_start) / 8;
state->sumBucketsOfX[bucket] += pixel;
state->sumBucketsOfY[bucket] += other_pixel;
state->sum2BucketsOfX[bucket] += pixel * pixel;
state->sum2BucketsOfY[bucket] += other_pixel * other_pixel;
state->sum2Buckets[bucket] += pixel * other_pixel;
}
c1 = (COLOR_Y_MAX * 0.01f) * (COLOR_Y_MAX * 0.01f);
c2 = (COLOR_Y_MAX * 0.03f) * (COLOR_Y_MAX * 0.03f);
break;
}
default: {
break;
}
}
// https://en.wikipedia.org/wiki/Structural_similarity
if ((!((state->lines_processed + 1) % 8)) || ((state->lines_processed + 1) == state->lines)) {
for (x = x_start; x < x_end; x += 8) {
int bucket = (x - x_start) / 8;
int w = IM_MIN((x_end - x), 8);
int h = IM_MIN((state->lines - state->lines_processed), 8);
float size = w * h;
// Dividng the sum squared buckets by size causes a loss of accuracy which results in
// the single pass standard deviation formula giving the wrong answer. To bypass this
// vx, vy, vxy have been multiplied by size which will be divided back out in the final
// ssim calculation (given c1/c2 ~= 0).
float mx = state->sumBucketsOfX[bucket] / size;
float my = state->sumBucketsOfY[bucket] / size;
float vx = state->sum2BucketsOfX[bucket] - (size * mx * mx);
float vy = state->sum2BucketsOfY[bucket] - (size * my * my);
float vxy = state->sum2Buckets[bucket] - (size * mx * my);
float ssim = (((2 * mx * my) + c1) * ((2 * vxy) + c2)) /
(((mx * mx) + (my * my) + c1) * (vx + vy + c2));
if (state->dssim) {
ssim = (1.0f - ssim) / 2.0f;
}
state->similarity_sum += ssim;
state->similarity_sum_2 += ssim * ssim;
state->similarity_min = IM_MIN(state->similarity_min, ssim);
state->similarity_max = IM_MAX(state->similarity_max, ssim);
state->sumBucketsOfX[bucket] = 0;
state->sumBucketsOfY[bucket] = 0;
state->sum2BucketsOfX[bucket] = 0;
state->sum2BucketsOfY[bucket] = 0;
state->sum2Buckets[bucket] = 0;
}
}
state->lines_processed += 1;
}
void imlib_get_similarity(image_t *img,
image_t *other,
int x_start,
int y_start,
float x_scale,
float y_scale,
rectangle_t *roi,
int rgb_channel,
int alpha,
const uint16_t *color_palette,
const uint8_t *alpha_palette,
image_hint_t hint,
bool dssim,
float *avg,
float *std,
float *min,
float *max) {
point_t p0, p1;
imlib_draw_image_get_bounds(img, other, x_start, y_start, x_scale, y_scale, roi,
alpha, alpha_palette, hint, &p0, &p1);
int h_blocks = ((p1.x - p0.x) + 7) / 8;
int v_blocks = ((p1.y - p0.y) + 7) / 8;
int blocks = h_blocks * v_blocks;
if (!blocks) {
return;
}
int int_h_blocks = h_blocks * sizeof(int);
imlib_similarity_line_op_state_t state;
state.dssim = dssim;
state.sumBucketsOfX = fb_alloc0(int_h_blocks * 5, FB_ALLOC_NO_HINT);
state.sumBucketsOfY = state.sumBucketsOfX + int_h_blocks;
state.sum2BucketsOfX = state.sumBucketsOfY + int_h_blocks;
state.sum2BucketsOfY = state.sum2BucketsOfX + int_h_blocks;
state.sum2Buckets = state.sum2BucketsOfY + int_h_blocks;
state.similarity_sum = 0.0f;
state.similarity_sum_2 = 0.0f;
state.similarity_min = FLT_MAX;
state.similarity_max = -FLT_MAX;
state.lines_processed = 0;
state.lines = p1.y - p0.y;
void *dst_row_override = fb_alloc0(image_line_size(img), FB_ALLOC_CACHE_ALIGN);
imlib_draw_image(img, other, x_start, y_start, x_scale, y_scale, roi,
rgb_channel, alpha, color_palette, alpha_palette, hint,
NULL, imlib_similarity_line_op, &state, dst_row_override);
*avg = state.similarity_sum / blocks;
*std = fast_sqrtf((state.similarity_sum_2 / blocks) - ((*avg) * (*avg)));
*min = state.similarity_min;
*max = state.similarity_max;
fb_free(); // dst_row_override
fb_free(); // sumBucketsOfX
}
#endif // IMLIB_ENABLE_GET_SIMILARITY
void imlib_get_histogram(histogram_t *out, image_t *ptr, rectangle_t *roi, list_t *thresholds, bool invert, image_t *other) {
switch (ptr->pixfmt) {
case PIXFORMAT_BINARY: {
memset(out->LBins, 0, out->LBinCount * sizeof(uint32_t));
int pixel_count = roi->w * roi->h;
float mult = (out->LBinCount - 1) / ((float) (COLOR_BINARY_MAX - COLOR_BINARY_MIN));
if ((!thresholds) || (!list_size(thresholds))) {
// Fast histogram code when no color thresholds list...
if (!other) {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x);
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_BINARY_MIN) * mult)]++;
}
}
} else {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x) ^ IMAGE_GET_BINARY_PIXEL_FAST(other_row_ptr, x);
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_BINARY_MIN) * mult)]++;
}
}
}
} else {
// Reset pixel count.
pixel_count = 0;
if (!other) {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x);
if (COLOR_THRESHOLD_BINARY(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_BINARY_MIN) * mult)]++;
pixel_count++;
}
}
}
}
} else {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x) ^ IMAGE_GET_BINARY_PIXEL_FAST(other_row_ptr,
x);
if (COLOR_THRESHOLD_BINARY(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_BINARY_MIN) * mult)]++;
pixel_count++;
}
}
}
}
}
}
float pixels = IM_DIV(1, ((float) pixel_count));
for (int i = 0, j = out->LBinCount; i < j; i++) {
out->LBins[i] = ((uint32_t *) out->LBins)[i] * pixels;
}
break;
}
case PIXFORMAT_GRAYSCALE: {
memset(out->LBins, 0, out->LBinCount * sizeof(uint32_t));
int pixel_count = roi->w * roi->h;
float mult = (out->LBinCount - 1) / ((float) (COLOR_GRAYSCALE_MAX - COLOR_GRAYSCALE_MIN));
if ((!thresholds) || (!list_size(thresholds))) {
// Fast histogram code when no color thresholds list...
if (!other) {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x);
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_GRAYSCALE_MIN) * mult)]++;
}
}
} else {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel =
abs(IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x) - IMAGE_GET_GRAYSCALE_PIXEL_FAST(other_row_ptr,
x));
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_GRAYSCALE_MIN) * mult)]++;
}
}
}
} else {
// Reset pixel count.
pixel_count = 0;
if (!other) {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x);
if (COLOR_THRESHOLD_GRAYSCALE(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_GRAYSCALE_MIN) * mult)]++;
pixel_count++;
}
}
}
}
} else {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel =
abs(IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr,
x) - IMAGE_GET_GRAYSCALE_PIXEL_FAST(other_row_ptr,
x));
if (COLOR_THRESHOLD_GRAYSCALE(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((pixel - COLOR_GRAYSCALE_MIN) * mult)]++;
pixel_count++;
}
}
}
}
}
}
float pixels = IM_DIV(1, ((float) pixel_count));
for (int i = 0, j = out->LBinCount; i < j; i++) {
out->LBins[i] = ((uint32_t *) out->LBins)[i] * pixels;
}
break;
}
case PIXFORMAT_RGB565: {
memset(out->LBins, 0, out->LBinCount * sizeof(uint32_t));
memset(out->ABins, 0, out->ABinCount * sizeof(uint32_t));
memset(out->BBins, 0, out->BBinCount * sizeof(uint32_t));
int pixel_count = roi->w * roi->h;
float l_mult = (out->LBinCount - 1) / ((float) (COLOR_L_MAX - COLOR_L_MIN));
float a_mult = (out->ABinCount - 1) / ((float) (COLOR_A_MAX - COLOR_A_MIN));
float b_mult = (out->BBinCount - 1) / ((float) (COLOR_B_MAX - COLOR_B_MIN));
if ((!thresholds) || (!list_size(thresholds))) {
// Fast histogram code when no color thresholds list...
if (!other) {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x);
((uint32_t *) out->LBins)[fast_roundf((COLOR_RGB565_TO_L(pixel) - COLOR_L_MIN) * l_mult)]++;
((uint32_t *) out->ABins)[fast_roundf((COLOR_RGB565_TO_A(pixel) - COLOR_A_MIN) * a_mult)]++;
((uint32_t *) out->BBins)[fast_roundf((COLOR_RGB565_TO_B(pixel) - COLOR_B_MIN) * b_mult)]++;
}
}
} else {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x);
int other_pixel = IMAGE_GET_RGB565_PIXEL_FAST(other_row_ptr, x);
int r = abs(COLOR_RGB565_TO_R5(pixel) - COLOR_RGB565_TO_R5(other_pixel));
int g = abs(COLOR_RGB565_TO_G6(pixel) - COLOR_RGB565_TO_G6(other_pixel));
int b = abs(COLOR_RGB565_TO_B5(pixel) - COLOR_RGB565_TO_B5(other_pixel));
pixel = COLOR_R5_G6_B5_TO_RGB565(r, g, b);
((uint32_t *) out->LBins)[fast_roundf((COLOR_RGB565_TO_L(pixel) - COLOR_L_MIN) * l_mult)]++;
((uint32_t *) out->ABins)[fast_roundf((COLOR_RGB565_TO_A(pixel) - COLOR_A_MIN) * a_mult)]++;
((uint32_t *) out->BBins)[fast_roundf((COLOR_RGB565_TO_B(pixel) - COLOR_B_MIN) * b_mult)]++;
}
}
}
} else {
// Reset pixel count.
pixel_count = 0;
if (!other) {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x);
if (COLOR_THRESHOLD_RGB565(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((COLOR_RGB565_TO_L(pixel) - COLOR_L_MIN) * l_mult)]++;
((uint32_t *) out->ABins)[fast_roundf((COLOR_RGB565_TO_A(pixel) - COLOR_A_MIN) * a_mult)]++;
((uint32_t *) out->BBins)[fast_roundf((COLOR_RGB565_TO_B(pixel) - COLOR_B_MIN) * b_mult)]++;
pixel_count++;
}
}
}
}
} else {
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y++) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y),
*other_row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(other, y);
for (int x = roi->x, xx = roi->x + roi->w; x < xx; x++) {
int pixel = IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x);
int other_pixel = IMAGE_GET_RGB565_PIXEL_FAST(other_row_ptr, x);
int r = abs(COLOR_RGB565_TO_R5(pixel) - COLOR_RGB565_TO_R5(other_pixel));
int g = abs(COLOR_RGB565_TO_G6(pixel) - COLOR_RGB565_TO_G6(other_pixel));
int b = abs(COLOR_RGB565_TO_B5(pixel) - COLOR_RGB565_TO_B5(other_pixel));
pixel = COLOR_R5_G6_B5_TO_RGB565(r, g, b);
if (COLOR_THRESHOLD_RGB565(pixel, lnk_data, invert)) {
((uint32_t *) out->LBins)[fast_roundf((COLOR_RGB565_TO_L(pixel) - COLOR_L_MIN) * l_mult)]++;
((uint32_t *) out->ABins)[fast_roundf((COLOR_RGB565_TO_A(pixel) - COLOR_A_MIN) * a_mult)]++;
((uint32_t *) out->BBins)[fast_roundf((COLOR_RGB565_TO_B(pixel) - COLOR_B_MIN) * b_mult)]++;
pixel_count++;
}
}
}
}
}
}
float pixels = IM_DIV(1, ((float) pixel_count));
for (int i = 0, j = out->LBinCount; i < j; i++) {
out->LBins[i] = ((uint32_t *) out->LBins)[i] * pixels;
}
for (int i = 0, j = out->ABinCount; i < j; i++) {
out->ABins[i] = ((uint32_t *) out->ABins)[i] * pixels;
}
for (int i = 0, j = out->BBinCount; i < j; i++) {
out->BBins[i] = ((uint32_t *) out->BBins)[i] * pixels;
}
break;
}
default: {
break;
}
}
}
void imlib_get_percentile(percentile_t *out, pixformat_t pixfmt, histogram_t *ptr, float percentile) {
memset(out, 0, sizeof(percentile_t));
switch (pixfmt) {
case PIXFORMAT_BINARY: {
float mult = (COLOR_BINARY_MAX - COLOR_BINARY_MIN) / ((float) (ptr->LBinCount - 1));
float median_count = 0;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
if ((median_count < percentile) && (percentile <= (median_count + ptr->LBins[i]))) {
out->LValue = fast_floorf((i * mult) + COLOR_BINARY_MIN);
break;
}
median_count += ptr->LBins[i];
}
break;
}
case PIXFORMAT_GRAYSCALE: {
float mult = (COLOR_GRAYSCALE_MAX - COLOR_GRAYSCALE_MIN) / ((float) (ptr->LBinCount - 1));
float median_count = 0;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
if ((median_count < percentile) && (percentile <= (median_count + ptr->LBins[i]))) {
out->LValue = fast_floorf((i * mult) + COLOR_GRAYSCALE_MIN);
break;
}
median_count += ptr->LBins[i];
}
break;
}
case PIXFORMAT_RGB565: {
{
float mult = (COLOR_L_MAX - COLOR_L_MIN) / ((float) (ptr->LBinCount - 1));
float median_count = 0;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
if ((median_count < percentile) && (percentile <= (median_count + ptr->LBins[i]))) {
out->LValue = fast_floorf((i * mult) + COLOR_L_MIN);
break;
}
median_count += ptr->LBins[i];
}
}
{
float mult = (COLOR_A_MAX - COLOR_A_MIN) / ((float) (ptr->ABinCount - 1));
float median_count = 0;
for (int i = 0, j = ptr->ABinCount; i < j; i++) {
if ((median_count < percentile) && (percentile <= (median_count + ptr->ABins[i]))) {
out->AValue = fast_floorf((i * mult) + COLOR_A_MIN);
break;
}
median_count += ptr->ABins[i];
}
}
{
float mult = (COLOR_B_MAX - COLOR_B_MIN) / ((float) (ptr->BBinCount - 1));
float median_count = 0;
for (int i = 0, j = ptr->BBinCount; i < j; i++) {
if ((median_count < percentile) && (percentile <= (median_count + ptr->BBins[i]))) {
out->BValue = fast_floorf((i * mult) + COLOR_B_MIN);
break;
}
median_count += ptr->BBins[i];
}
}
break;
}
default: {
break;
}
}
}
static int ostu(int bincount, float *bins) {
float cdf[bincount]; memset(cdf, 0, bincount * sizeof(float));
float weighted_cdf[bincount]; memset(weighted_cdf, 0, bincount * sizeof(float));
cdf[0] = bins[0];
weighted_cdf[0] = 0 * bins[0];
for (int i = 1; i < bincount; i++) {
cdf[i] = cdf[i - 1] + bins[i];
weighted_cdf[i] = weighted_cdf[i - 1] + (i * bins[i]);
}
float variance[bincount]; memset(variance, 0, bincount * sizeof(float));
float max_variance = 0.0f;
int threshold = 0;
for (int i = 0, ii = bincount - 1; i < ii; i++) {
if ((cdf[i] != 0.0f) && (cdf[i] != 1.0f)) {
variance[i] = powf((cdf[i] * weighted_cdf[bincount - 1]) - weighted_cdf[i], 2.0f) / (cdf[i] * (1.0f - cdf[i]));
} else {
variance[i] = 0.0f;
}
if (variance[i] > max_variance) {
max_variance = variance[i];
threshold = i;
}
}
return threshold;
}
void imlib_get_threshold(threshold_t *out, pixformat_t pixfmt, histogram_t *ptr) {
memset(out, 0, sizeof(threshold_t));
switch (pixfmt) {
case PIXFORMAT_BINARY: {
out->LValue = (ostu(ptr->LBinCount, ptr->LBins) * (COLOR_BINARY_MAX - COLOR_BINARY_MIN)) / (ptr->LBinCount - 1);
break;
}
case PIXFORMAT_GRAYSCALE: {
out->LValue =
(ostu(ptr->LBinCount, ptr->LBins) * (COLOR_GRAYSCALE_MAX - COLOR_GRAYSCALE_MIN)) / (ptr->LBinCount - 1);
break;
}
case PIXFORMAT_RGB565: {
out->LValue = (ostu(ptr->LBinCount, ptr->LBins) * (COLOR_L_MAX - COLOR_L_MIN)) / (ptr->LBinCount - 1);
out->AValue = ((ostu(ptr->ABinCount, ptr->ABins) * (COLOR_A_MAX - COLOR_A_MIN)) / (ptr->ABinCount - 1)) +
COLOR_A_MIN;
out->BValue = ((ostu(ptr->BBinCount, ptr->BBins) * (COLOR_B_MAX - COLOR_B_MIN)) / (ptr->BBinCount - 1)) +
COLOR_B_MIN;
break;
}
default: {
break;
}
}
}
void imlib_get_statistics(statistics_t *out, pixformat_t pixfmt, histogram_t *ptr) {
memset(out, 0, sizeof(statistics_t));
switch (pixfmt) {
case PIXFORMAT_BINARY: {
float mult = (COLOR_BINARY_MAX - COLOR_BINARY_MIN) / ((float) (ptr->LBinCount - 1));
float avg = 0;
float stdev = 0;
float median_count = 0;
float mode_count = 0;
bool min_flag = false;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
float value_f = (i * mult) + COLOR_BINARY_MIN;
int value = fast_floorf(value_f);
avg += value_f * ptr->LBins[i];
stdev += value_f * value_f * ptr->LBins[i];
if ((median_count < 0.25f) && (0.25f <= (median_count + ptr->LBins[i]))) {
out->LLQ = value;
}
if ((median_count < 0.5f) && (0.5f <= (median_count + ptr->LBins[i]))) {
out->LMedian = value;
}
if ((median_count < 0.75f) && (0.75f <= (median_count + ptr->LBins[i]))) {
out->LUQ = value;
}
if (ptr->LBins[i] > mode_count) {
mode_count = ptr->LBins[i];
out->LMode = value;
}
if ((ptr->LBins[i] > 0.0f) && (!min_flag)) {
min_flag = true;
out->LMin = value;
}
if (ptr->LBins[i] > 0.0f) {
out->LMax = value;
}
median_count += ptr->LBins[i];
}
out->LMean = fast_floorf(avg);
out->LSTDev = fast_floorf(fast_sqrtf(stdev - (avg * avg)));
break;
}
case PIXFORMAT_GRAYSCALE: {
float mult = (COLOR_GRAYSCALE_MAX - COLOR_GRAYSCALE_MIN) / ((float) (ptr->LBinCount - 1));
float avg = 0;
float stdev = 0;
float median_count = 0;
float mode_count = 0;
bool min_flag = false;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
float value_f = (i * mult) + COLOR_GRAYSCALE_MIN;
int value = fast_floorf(value_f);
avg += value_f * ptr->LBins[i];
stdev += value_f * value_f * ptr->LBins[i];
if ((median_count < 0.25f) && (0.25f <= (median_count + ptr->LBins[i]))) {
out->LLQ = value;
}
if ((median_count < 0.5f) && (0.5f <= (median_count + ptr->LBins[i]))) {
out->LMedian = value;
}
if ((median_count < 0.75f) && (0.75f <= (median_count + ptr->LBins[i]))) {
out->LUQ = value;
}
if (ptr->LBins[i] > mode_count) {
mode_count = ptr->LBins[i];
out->LMode = value;
}
if ((ptr->LBins[i] > 0.0f) && (!min_flag)) {
min_flag = true;
out->LMin = value;
}
if (ptr->LBins[i] > 0.0f) {
out->LMax = value;
}
median_count += ptr->LBins[i];
}
out->LMean = fast_floorf(avg);
out->LSTDev = fast_floorf(fast_sqrtf(stdev - (avg * avg)));
break;
}
case PIXFORMAT_RGB565: {
{
float mult = (COLOR_L_MAX - COLOR_L_MIN) / ((float) (ptr->LBinCount - 1));
float avg = 0;
float stdev = 0;
float median_count = 0;
float mode_count = 0;
bool min_flag = false;
for (int i = 0, j = ptr->LBinCount; i < j; i++) {
float value_f = (i * mult) + COLOR_L_MIN;
int value = fast_floorf(value_f);
avg += value_f * ptr->LBins[i];
stdev += value_f * value_f * ptr->LBins[i];
if ((median_count < 0.25f) && (0.25f <= (median_count + ptr->LBins[i]))) {
out->LLQ = value;
}
if ((median_count < 0.5f) && (0.5f <= (median_count + ptr->LBins[i]))) {
out->LMedian = value;
}
if ((median_count < 0.75f) && (0.75f <= (median_count + ptr->LBins[i]))) {
out->LUQ = value;
}
if (ptr->LBins[i] > mode_count) {
mode_count = ptr->LBins[i];
out->LMode = value;
}
if ((ptr->LBins[i] > 0.0f) && (!min_flag)) {
min_flag = true;
out->LMin = value;
}
if (ptr->LBins[i] > 0.0f) {
out->LMax = value;
}
median_count += ptr->LBins[i];
}
out->LMean = fast_floorf(avg);
out->LSTDev = fast_floorf(fast_sqrtf(stdev - (avg * avg)));
}
{
float mult = (COLOR_A_MAX - COLOR_A_MIN) / ((float) (ptr->ABinCount - 1));
float avg = 0;
float stdev = 0;
float median_count = 0;
float mode_count = 0;
bool min_flag = false;
for (int i = 0, j = ptr->ABinCount; i < j; i++) {
float value_f = (i * mult) + COLOR_A_MIN;
int value = fast_floorf(value_f);
avg += value_f * ptr->ABins[i];
stdev += value_f * value_f * ptr->ABins[i];
if ((median_count < 0.25f) && (0.25f <= (median_count + ptr->ABins[i]))) {
out->ALQ = value;
}
if ((median_count < 0.5f) && (0.5f <= (median_count + ptr->ABins[i]))) {
out->AMedian = value;
}
if ((median_count < 0.75f) && (0.75f <= (median_count + ptr->ABins[i]))) {
out->AUQ = value;
}
if (ptr->ABins[i] > mode_count) {
mode_count = ptr->ABins[i];
out->AMode = value;
}
if ((ptr->ABins[i] > 0.0f) && (!min_flag)) {
min_flag = true;
out->AMin = value;
}
if (ptr->ABins[i] > 0.0f) {
out->AMax = value;
}
median_count += ptr->ABins[i];
}
out->AMean = fast_floorf(avg);
out->ASTDev = fast_floorf(fast_sqrtf(stdev - (avg * avg)));
}
{
float mult = (COLOR_B_MAX - COLOR_B_MIN) / ((float) (ptr->BBinCount - 1));
float avg = 0;
float stdev = 0;
float median_count = 0;
float mode_count = 0;
bool min_flag = false;
for (int i = 0, j = ptr->BBinCount; i < j; i++) {
float value_f = (i * mult) + COLOR_B_MIN;
int value = fast_floorf(value_f);
avg += value_f * ptr->BBins[i];
stdev += value_f * value_f * ptr->BBins[i];
if ((median_count < 0.25f) && (0.25f <= (median_count + ptr->BBins[i]))) {
out->BLQ = value;
}
if ((median_count < 0.5f) && (0.5f <= (median_count + ptr->BBins[i]))) {
out->BMedian = value;
}
if ((median_count < 0.75f) && (0.75f <= (median_count + ptr->BBins[i]))) {
out->BUQ = value;
}
if (ptr->BBins[i] > mode_count) {
mode_count = ptr->BBins[i];
out->BMode = value;
}
if ((ptr->BBins[i] > 0.0f) && (!min_flag)) {
min_flag = true;
out->BMin = value;
}
if (ptr->BBins[i] > 0.0f) {
out->BMax = value;
}
median_count += ptr->BBins[i];
}
out->BMean = fast_floorf(avg);
out->BSTDev = fast_floorf(fast_sqrtf(stdev - (avg * avg)));
}
break;
}
default: {
break;
}
}
}
static int get_median(int *array, int array_sum, int array_len) {
const int median_threshold = (array_sum + 1) / 2;
int median_count = 0;
for (int i = 0; i < array_len; i++) {
if ((median_count < median_threshold) && (median_threshold <= (median_count + array[i]))) {
return i;
}
median_count += array[i];
}
return array_len - 1;
}
static int get_median_l(long long *array, long long array_sum, int array_len) {
const long long median_threshold = (array_sum + 1) / 2;
long long median_count = 0;
for (int i = 0; i < array_len; i++) {
if ((median_count < median_threshold) && (median_threshold <= (median_count + array[i]))) {
return i;
}
median_count += array[i];
}
return array_len - 1;
}
bool imlib_get_regression(find_lines_list_lnk_data_t *out,
image_t *ptr,
rectangle_t *roi,
unsigned int x_stride,
unsigned int y_stride,
list_t *thresholds,
bool invert,
unsigned int area_threshold,
unsigned int pixels_threshold,
bool robust) {
bool result = false;
memset(out, 0, sizeof(find_lines_list_lnk_data_t));
if (!robust) {
// Least Squares
int blob_x1 = roi->x + roi->w - 1;
int blob_y1 = roi->y + roi->h - 1;
int blob_x2 = roi->x;
int blob_y2 = roi->y;
int blob_pixels = 0;
int blob_cx = 0;
int blob_cy = 0;
long long blob_a = 0;
long long blob_b = 0;
long long blob_c = 0;
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
switch (ptr->pixfmt) {
case PIXFORMAT_BINARY: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_BINARY(IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
blob_cx += x;
blob_cy += y;
blob_a += x * x;
blob_b += x * y;
blob_c += y * y;
}
}
}
break;
}
case PIXFORMAT_GRAYSCALE: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_GRAYSCALE(IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
blob_cx += x;
blob_cy += y;
blob_a += x * x;
blob_b += x * y;
blob_c += y * y;
}
}
}
break;
}
case PIXFORMAT_RGB565: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_RGB565(IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
blob_cx += x;
blob_cy += y;
blob_a += x * x;
blob_b += x * y;
blob_c += y * y;
}
}
}
break;
}
default: {
break;
}
}
}
int w = blob_x2 - blob_x1;
int h = blob_y2 - blob_y1;
if (blob_pixels && ((w * h) >= area_threshold) && (blob_pixels >= pixels_threshold)) {
// http://www.cse.usf.edu/~r1k/MachineVisionBook/MachineVision.files/MachineVision_Chapter2.pdf
// https://www.strchr.com/standard_deviation_in_one_pass
//
// a = sigma(x*x) + (mx*sigma(x)) + (mx*sigma(x)) + (sigma()*mx*mx)
// b = sigma(x*y) + (mx*sigma(y)) + (my*sigma(x)) + (sigma()*mx*my)
// c = sigma(y*y) + (my*sigma(y)) + (my*sigma(y)) + (sigma()*my*my)
//
// blob_a = sigma(x*x)
// blob_b = sigma(x*y)
// blob_c = sigma(y*y)
// blob_cx = sigma(x)
// blob_cy = sigma(y)
// blob_pixels = sigma()
int mx = blob_cx / blob_pixels; // x centroid
int my = blob_cy / blob_pixels; // y centroid
int small_blob_a = blob_a - ((mx * blob_cx) + (mx * blob_cx)) + (blob_pixels * mx * mx);
int small_blob_b = blob_b - ((mx * blob_cy) + (my * blob_cx)) + (blob_pixels * mx * my);
int small_blob_c = blob_c - ((my * blob_cy) + (my * blob_cy)) + (blob_pixels * my * my);
float rotation =
((small_blob_a !=
small_blob_c) ? (fast_atan2f(2 * small_blob_b, small_blob_a - small_blob_c) / 2.0f) : 1.570796f) + 1.570796f; // PI/2
out->theta = fast_roundf(rotation * 57.295780) % 180; // * (180 / PI)
if (out->theta < 0) {
out->theta += 180;
}
out->rho = fast_roundf(((mx - roi->x) * cos_table[out->theta]) + ((my - roi->y) * sin_table[out->theta]));
float part0 = (small_blob_a + small_blob_c) / 2.0f;
float f_b = (float) small_blob_b;
float f_a_c = (float) (small_blob_a - small_blob_c);
float part1 = fast_sqrtf((4 * f_b * f_b) + (f_a_c * f_a_c)) / 2.0f;
float p_add = fast_sqrtf(part0 + part1);
float p_sub = fast_sqrtf(part0 - part1);
float e_min = IM_MIN(p_add, p_sub);
float e_max = IM_MAX(p_add, p_sub);
out->magnitude = fast_roundf(e_max / e_min) - 1; // Circle -> [0, INF) -> Line
if ((45 <= out->theta) && (out->theta < 135)) {
// y = (r - x cos(t)) / sin(t)
out->line.x1 = 0;
out->line.y1 = fast_roundf((out->rho - (out->line.x1 * cos_table[out->theta])) / sin_table[out->theta]);
out->line.x2 = roi->w - 1;
out->line.y2 = fast_roundf((out->rho - (out->line.x2 * cos_table[out->theta])) / sin_table[out->theta]);
} else {
// x = (r - y sin(t)) / cos(t);
out->line.y1 = 0;
out->line.x1 = fast_roundf((out->rho - (out->line.y1 * sin_table[out->theta])) / cos_table[out->theta]);
out->line.y2 = roi->h - 1;
out->line.x2 = fast_roundf((out->rho - (out->line.y2 * sin_table[out->theta])) / cos_table[out->theta]);
}
if (lb_clip_line(&out->line, 0, 0, roi->w, roi->h)) {
out->line.x1 += roi->x;
out->line.y1 += roi->y;
out->line.x2 += roi->x;
out->line.y2 += roi->y;
// Move rho too.
out->rho += fast_roundf((roi->x * cos_table[out->theta]) + (roi->y * sin_table[out->theta]));
result = true;
} else {
memset(out, 0, sizeof(find_lines_list_lnk_data_t));
}
}
} else {
// Theil-Sen Estimator
int *x_histogram = fb_alloc0(ptr->w * sizeof(int), FB_ALLOC_NO_HINT);
int *y_histogram = fb_alloc0(ptr->h * sizeof(int), FB_ALLOC_NO_HINT);
long long *x_delta_histogram = fb_alloc0((2 * ptr->w) * sizeof(long long), FB_ALLOC_NO_HINT);
long long *y_delta_histogram = fb_alloc0((2 * ptr->h) * sizeof(long long), FB_ALLOC_NO_HINT);
uint32_t size;
point_t *points = (point_t *) fb_alloc_all(&size, FB_ALLOC_NO_HINT);
size_t points_max = size / sizeof(point_t);
size_t points_count = 0;
if (points_max) {
int blob_x1 = roi->x + roi->w - 1;
int blob_y1 = roi->y + roi->h - 1;
int blob_x2 = roi->x;
int blob_y2 = roi->y;
int blob_pixels = 0;
list_for_each(it, thresholds) {
color_thresholds_list_lnk_data_t *lnk_data = list_get_data(it);
switch (ptr->pixfmt) {
case PIXFORMAT_BINARY: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint32_t *row_ptr = IMAGE_COMPUTE_BINARY_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_BINARY(IMAGE_GET_BINARY_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
x_histogram[x]++;
y_histogram[y]++;
if (points_count < points_max) {
point_init(&points[points_count], x, y);
points_count += 1;
}
}
}
}
break;
}
case PIXFORMAT_GRAYSCALE: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint8_t *row_ptr = IMAGE_COMPUTE_GRAYSCALE_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_GRAYSCALE(IMAGE_GET_GRAYSCALE_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
x_histogram[x]++;
y_histogram[y]++;
if (points_count < points_max) {
point_init(&points[points_count], x, y);
points_count += 1;
}
}
}
}
break;
}
case PIXFORMAT_RGB565: {
for (int y = roi->y, yy = roi->y + roi->h; y < yy; y += y_stride) {
uint16_t *row_ptr = IMAGE_COMPUTE_RGB565_PIXEL_ROW_PTR(ptr, y);
for (int x = roi->x + (y % x_stride), xx = roi->x + roi->w; x < xx; x += x_stride) {
if (COLOR_THRESHOLD_RGB565(IMAGE_GET_RGB565_PIXEL_FAST(row_ptr, x), lnk_data, invert)) {
blob_x1 = IM_MIN(blob_x1, x);
blob_y1 = IM_MIN(blob_y1, y);
blob_x2 = IM_MAX(blob_x2, x);
blob_y2 = IM_MAX(blob_y2, y);
blob_pixels += 1;
x_histogram[x]++;
y_histogram[y]++;
if (points_count < points_max) {
point_init(&points[points_count], x, y);
points_count += 1;
}
}
}
}
break;
}
default: {
break;
}
}
}
int w = blob_x2 - blob_x1;
int h = blob_y2 - blob_y1;
if (blob_pixels && ((w * h) >= area_threshold) && (blob_pixels >= pixels_threshold)) {
long long delta_sum = (points_count * (points_count - 1)) / 2;
if (delta_sum) {
// The code below computes the average slope between all pairs of points.
// This is a N^2 operation that can easily blow up if the image is not threshold carefully...
for (int i = 0; i < points_count; i++) {
point_t *p0 = &points[i];
for (int j = i + 1; j < points_count; j++) {
point_t *p1 = &points[j];
// Note we allocated 1 extra above so we can do ptr->w instead of (ptr->w-1).
x_delta_histogram[p0->x - p1->x + ptr->w]++;
// Note we allocated 1 extra above so we can do ptr->h instead of (ptr->h-1).
y_delta_histogram[p0->y - p1->y + ptr->h]++;
}
}
int mx = get_median(x_histogram, blob_pixels, ptr->w); // Output doesn't need adjustment.
int my = get_median(y_histogram, blob_pixels, ptr->h); // Output doesn't need adjustment.
int mdx = get_median_l(x_delta_histogram, delta_sum, 2 * ptr->w) - ptr->w; // Fix offset.
int mdy = get_median_l(y_delta_histogram, delta_sum, 2 * ptr->h) - ptr->h; // Fix offset.
float rotation = (mdx ? fast_atan2f(mdy, mdx) : 1.570796f) + 1.570796f; // PI/2
out->theta = fast_roundf(rotation * 57.295780) % 180; // * (180 / PI)
if (out->theta < 0) {
out->theta += 180;
}
out->rho = fast_roundf(((mx - roi->x) * cos_table[out->theta]) + ((my - roi->y) * sin_table[out->theta]));
out->magnitude = fast_roundf(fast_sqrtf((mdx * mdx) + (mdy * mdy)));
if ((45 <= out->theta) && (out->theta < 135)) {
// y = (r - x cos(t)) / sin(t)
out->line.x1 = 0;
out->line.y1 = fast_roundf((out->rho - (out->line.x1 * cos_table[out->theta])) / sin_table[out->theta]);
out->line.x2 = roi->w - 1;
out->line.y2 = fast_roundf((out->rho - (out->line.x2 * cos_table[out->theta])) / sin_table[out->theta]);
} else {
// x = (r - y sin(t)) / cos(t);
out->line.y1 = 0;
out->line.x1 = fast_roundf((out->rho - (out->line.y1 * sin_table[out->theta])) / cos_table[out->theta]);
out->line.y2 = roi->h - 1;
out->line.x2 = fast_roundf((out->rho - (out->line.y2 * sin_table[out->theta])) / cos_table[out->theta]);
}
if (lb_clip_line(&out->line, 0, 0, roi->w, roi->h)) {
out->line.x1 += roi->x;
out->line.y1 += roi->y;
out->line.x2 += roi->x;
out->line.y2 += roi->y;
// Move rho too.
out->rho += fast_roundf((roi->x * cos_table[out->theta]) + (roi->y * sin_table[out->theta]));
result = true;
} else {
memset(out, 0, sizeof(find_lines_list_lnk_data_t));
}
}
}
}
fb_free(); // points
fb_free(); // y_delta_histogram
fb_free(); // x_delta_histogram
fb_free(); // y_histogram
fb_free(); // x_histogram
}
return result;
}