enc_detect_dots.cc (22534B)
1 // Copyright (c) the JPEG XL Project Authors. All rights reserved. 2 // 3 // Use of this source code is governed by a BSD-style 4 // license that can be found in the LICENSE file. 5 6 #include "lib/jxl/enc_detect_dots.h" 7 8 #include <stdint.h> 9 10 #include <algorithm> 11 #include <array> 12 #include <cmath> 13 #include <cstdio> 14 #include <utility> 15 #include <vector> 16 17 #undef HWY_TARGET_INCLUDE 18 #define HWY_TARGET_INCLUDE "lib/jxl/enc_detect_dots.cc" 19 #include <hwy/foreach_target.h> 20 #include <hwy/highway.h> 21 22 #include "lib/jxl/base/common.h" 23 #include "lib/jxl/base/compiler_specific.h" 24 #include "lib/jxl/base/data_parallel.h" 25 #include "lib/jxl/base/printf_macros.h" 26 #include "lib/jxl/base/status.h" 27 #include "lib/jxl/convolve.h" 28 #include "lib/jxl/enc_linalg.h" 29 #include "lib/jxl/image.h" 30 #include "lib/jxl/image_ops.h" 31 32 // Set JXL_DEBUG_DOT_DETECT to 1 to enable debugging. 33 #ifndef JXL_DEBUG_DOT_DETECT 34 #define JXL_DEBUG_DOT_DETECT 0 35 #endif 36 37 HWY_BEFORE_NAMESPACE(); 38 namespace jxl { 39 namespace HWY_NAMESPACE { 40 41 // These templates are not found via ADL. 42 using hwy::HWY_NAMESPACE::Add; 43 using hwy::HWY_NAMESPACE::Mul; 44 using hwy::HWY_NAMESPACE::Sub; 45 46 StatusOr<ImageF> SumOfSquareDifferences(const Image3F& forig, 47 const Image3F& smooth, 48 ThreadPool* pool) { 49 const HWY_FULL(float) d; 50 const auto color_coef0 = Set(d, 0.0f); 51 const auto color_coef1 = Set(d, 10.0f); 52 const auto color_coef2 = Set(d, 0.0f); 53 54 JXL_ASSIGN_OR_RETURN(ImageF sum_of_squares, 55 ImageF::Create(forig.xsize(), forig.ysize())); 56 JXL_CHECK(RunOnPool( 57 pool, 0, forig.ysize(), ThreadPool::NoInit, 58 [&](const uint32_t task, size_t thread) { 59 const size_t y = static_cast<size_t>(task); 60 const float* JXL_RESTRICT orig_row0 = forig.Plane(0).ConstRow(y); 61 const float* JXL_RESTRICT orig_row1 = forig.Plane(1).ConstRow(y); 62 const float* JXL_RESTRICT orig_row2 = forig.Plane(2).ConstRow(y); 63 const float* JXL_RESTRICT smooth_row0 = smooth.Plane(0).ConstRow(y); 64 const float* JXL_RESTRICT smooth_row1 = smooth.Plane(1).ConstRow(y); 65 const float* JXL_RESTRICT smooth_row2 = smooth.Plane(2).ConstRow(y); 66 float* JXL_RESTRICT sos_row = sum_of_squares.Row(y); 67 68 for (size_t x = 0; x < forig.xsize(); x += Lanes(d)) { 69 auto v0 = Sub(Load(d, orig_row0 + x), Load(d, smooth_row0 + x)); 70 auto v1 = Sub(Load(d, orig_row1 + x), Load(d, smooth_row1 + x)); 71 auto v2 = Sub(Load(d, orig_row2 + x), Load(d, smooth_row2 + x)); 72 v0 = Mul(Mul(v0, v0), color_coef0); 73 v1 = Mul(Mul(v1, v1), color_coef1); 74 v2 = Mul(Mul(v2, v2), color_coef2); 75 const auto sos = 76 Add(v0, Add(v1, v2)); // weighted sum of square diffs 77 Store(sos, d, sos_row + x); 78 } 79 }, 80 "ComputeEnergyImage")); 81 return sum_of_squares; 82 } 83 84 // NOLINTNEXTLINE(google-readability-namespace-comments) 85 } // namespace HWY_NAMESPACE 86 } // namespace jxl 87 HWY_AFTER_NAMESPACE(); 88 89 #if HWY_ONCE 90 namespace jxl { 91 HWY_EXPORT(SumOfSquareDifferences); // Local function 92 93 const int kEllipseWindowSize = 5; 94 95 namespace { 96 struct GaussianEllipse { 97 double x; // position in x 98 double y; // position in y 99 double sigma_x; // scale in x 100 double sigma_y; // scale in y 101 double angle; // ellipse rotation in radians 102 std::array<double, 3> intensity; // intensity in each channel 103 104 // The following variables do not need to be encoded 105 double l2_loss; // error after the Gaussian was fit 106 double l1_loss; 107 double ridge_loss; // the l2_loss plus regularization term 108 double custom_loss; // experimental custom loss 109 std::array<double, 3> bgColor; // best background color 110 size_t neg_pixels; // number of negative pixels when subtracting dot 111 std::array<double, 3> neg_value; // debt due to channel truncation 112 }; 113 double DotGaussianModel(double dx, double dy, double ct, double st, 114 double sigma_x, double sigma_y, double intensity) { 115 double rx = ct * dx + st * dy; 116 double ry = -st * dx + ct * dy; 117 double md = (rx * rx / sigma_x) + (ry * ry / sigma_y); 118 double value = intensity * exp(-0.5 * md); 119 return value; 120 } 121 122 constexpr bool kOptimizeBackground = true; 123 124 // Gaussian that smooths noise but preserves dots 125 const WeightsSeparable5& WeightsSeparable5Gaussian0_65() { 126 constexpr float w0 = 0.558311f; 127 constexpr float w1 = 0.210395f; 128 constexpr float w2 = 0.010449f; 129 static constexpr WeightsSeparable5 weights = { 130 {HWY_REP4(w0), HWY_REP4(w1), HWY_REP4(w2)}, 131 {HWY_REP4(w0), HWY_REP4(w1), HWY_REP4(w2)}}; 132 return weights; 133 } 134 135 // (Iterated) Gaussian that removes dots. 136 const WeightsSeparable5& WeightsSeparable5Gaussian3() { 137 constexpr float w0 = 0.222338f; 138 constexpr float w1 = 0.210431f; 139 constexpr float w2 = 0.1784f; 140 static constexpr WeightsSeparable5 weights = { 141 {HWY_REP4(w0), HWY_REP4(w1), HWY_REP4(w2)}, 142 {HWY_REP4(w0), HWY_REP4(w1), HWY_REP4(w2)}}; 143 return weights; 144 } 145 146 StatusOr<ImageF> ComputeEnergyImage(const Image3F& orig, Image3F* smooth, 147 ThreadPool* pool) { 148 // Prepare guidance images for dot selection. 149 JXL_ASSIGN_OR_RETURN(Image3F forig, 150 Image3F::Create(orig.xsize(), orig.ysize())); 151 JXL_ASSIGN_OR_RETURN(*smooth, Image3F::Create(orig.xsize(), orig.ysize())); 152 Rect rect(orig); 153 154 const auto& weights1 = WeightsSeparable5Gaussian0_65(); 155 const auto& weights3 = WeightsSeparable5Gaussian3(); 156 157 for (size_t c = 0; c < 3; ++c) { 158 // Use forig as temporary storage to reduce memory and keep it warmer. 159 Separable5(orig.Plane(c), rect, weights3, pool, &forig.Plane(c)); 160 Separable5(forig.Plane(c), rect, weights3, pool, &smooth->Plane(c)); 161 Separable5(orig.Plane(c), rect, weights1, pool, &forig.Plane(c)); 162 } 163 164 return HWY_DYNAMIC_DISPATCH(SumOfSquareDifferences)(forig, *smooth, pool); 165 } 166 167 struct Pixel { 168 int x; 169 int y; 170 }; 171 172 Pixel operator+(const Pixel& a, const Pixel& b) { 173 return Pixel{a.x + b.x, a.y + b.y}; 174 } 175 176 // Maximum area in pixels of a ellipse 177 const size_t kMaxCCSize = 1000; 178 179 // Extracts a connected component from a Binary image where seed is part 180 // of the component 181 bool ExtractComponent(const Rect& rect, ImageF* img, std::vector<Pixel>* pixels, 182 const Pixel& seed, double threshold) { 183 static const std::vector<Pixel> neighbors{{1, -1}, {1, 0}, {1, 1}, {0, -1}, 184 {0, 1}, {-1, -1}, {-1, 1}, {1, 0}}; 185 std::vector<Pixel> q{seed}; 186 while (!q.empty()) { 187 Pixel current = q.back(); 188 q.pop_back(); 189 pixels->push_back(current); 190 if (pixels->size() > kMaxCCSize) return false; 191 for (const Pixel& delta : neighbors) { 192 Pixel child = current + delta; 193 if (child.x >= 0 && static_cast<size_t>(child.x) < rect.xsize() && 194 child.y >= 0 && static_cast<size_t>(child.y) < rect.ysize()) { 195 float* value = &rect.Row(img, child.y)[child.x]; 196 if (*value > threshold) { 197 *value = 0.0; 198 q.push_back(child); 199 } 200 } 201 } 202 } 203 return true; 204 } 205 206 inline bool PointInRect(const Rect& r, const Pixel& p) { 207 return (static_cast<size_t>(p.x) >= r.x0() && 208 static_cast<size_t>(p.x) < (r.x0() + r.xsize()) && 209 static_cast<size_t>(p.y) >= r.y0() && 210 static_cast<size_t>(p.y) < (r.y0() + r.ysize())); 211 } 212 213 struct ConnectedComponent { 214 ConnectedComponent(const Rect& bounds, const std::vector<Pixel>&& pixels) 215 : bounds(bounds), pixels(pixels) {} 216 Rect bounds; 217 std::vector<Pixel> pixels; 218 float maxEnergy; 219 float meanEnergy; 220 float varEnergy; 221 float meanBg; 222 float varBg; 223 float score; 224 Pixel mode; 225 226 void CompStats(const ImageF& energy, const Rect& rect, int extra) { 227 maxEnergy = 0.0; 228 meanEnergy = 0.0; 229 varEnergy = 0.0; 230 meanBg = 0.0; 231 varBg = 0.0; 232 int nIn = 0; 233 int nOut = 0; 234 mode.x = 0; 235 mode.y = 0; 236 for (int sy = -extra; sy < (static_cast<int>(bounds.ysize()) + extra); 237 sy++) { 238 int y = sy + static_cast<int>(bounds.y0()); 239 if (y < 0 || static_cast<size_t>(y) >= rect.ysize()) continue; 240 const float* JXL_RESTRICT erow = rect.ConstRow(energy, y); 241 for (int sx = -extra; sx < (static_cast<int>(bounds.xsize()) + extra); 242 sx++) { 243 int x = sx + static_cast<int>(bounds.x0()); 244 if (x < 0 || static_cast<size_t>(x) >= rect.xsize()) continue; 245 if (erow[x] > maxEnergy) { 246 maxEnergy = erow[x]; 247 mode.x = x; 248 mode.y = y; 249 } 250 if (PointInRect(bounds, Pixel{x, y})) { 251 meanEnergy += erow[x]; 252 varEnergy += erow[x] * erow[x]; 253 nIn++; 254 } else { 255 meanBg += erow[x]; 256 varBg += erow[x] * erow[x]; 257 nOut++; 258 } 259 } 260 } 261 meanEnergy = meanEnergy / nIn; 262 meanBg = meanBg / nOut; 263 varEnergy = (varEnergy / nIn) - meanEnergy * meanEnergy; 264 varBg = (varBg / nOut) - meanBg * meanBg; 265 score = (meanEnergy - meanBg) / std::sqrt(varBg); 266 } 267 }; 268 269 Rect BoundingRectangle(const std::vector<Pixel>& pixels) { 270 JXL_ASSERT(!pixels.empty()); 271 int low_x; 272 int high_x; 273 int low_y; 274 int high_y; 275 low_x = high_x = pixels[0].x; 276 low_y = high_y = pixels[0].y; 277 for (const Pixel& p : pixels) { 278 low_x = std::min(low_x, p.x); 279 high_x = std::max(high_x, p.x); 280 low_y = std::min(low_y, p.y); 281 high_y = std::max(high_y, p.y); 282 } 283 return Rect(low_x, low_y, high_x - low_x + 1, high_y - low_y + 1); 284 } 285 286 StatusOr<std::vector<ConnectedComponent>> FindCC(const ImageF& energy, 287 const Rect& rect, double t_low, 288 double t_high, 289 uint32_t maxWindow, 290 double minScore) { 291 const int kExtraRect = 4; 292 JXL_ASSIGN_OR_RETURN(ImageF img, 293 ImageF::Create(energy.xsize(), energy.ysize())); 294 CopyImageTo(energy, &img); 295 std::vector<ConnectedComponent> ans; 296 for (size_t y = 0; y < rect.ysize(); y++) { 297 float* JXL_RESTRICT row = rect.Row(&img, y); 298 for (size_t x = 0; x < rect.xsize(); x++) { 299 if (row[x] > t_high) { 300 std::vector<Pixel> pixels; 301 row[x] = 0.0; 302 bool success = ExtractComponent( 303 rect, &img, &pixels, 304 Pixel{static_cast<int>(x), static_cast<int>(y)}, t_low); 305 if (!success) continue; 306 #if JXL_DEBUG_DOT_DETECT 307 for (size_t i = 0; i < pixels.size(); i++) { 308 fprintf(stderr, "(%d,%d) ", pixels[i].x, pixels[i].y); 309 } 310 fprintf(stderr, "\n"); 311 #endif // JXL_DEBUG_DOT_DETECT 312 Rect bounds = BoundingRectangle(pixels); 313 if (bounds.xsize() < maxWindow && bounds.ysize() < maxWindow) { 314 ConnectedComponent cc{bounds, std::move(pixels)}; 315 cc.CompStats(energy, rect, kExtraRect); 316 if (cc.score < minScore) continue; 317 JXL_DEBUG(JXL_DEBUG_DOT_DETECT, 318 "cc mode: (%d,%d), max: %f, bgMean: %f bgVar: " 319 "%f bound:(%" PRIuS ",%" PRIuS ",%" PRIuS ",%" PRIuS ")\n", 320 cc.mode.x, cc.mode.y, cc.maxEnergy, cc.meanEnergy, 321 cc.varEnergy, cc.bounds.x0(), cc.bounds.y0(), 322 cc.bounds.xsize(), cc.bounds.ysize()); 323 ans.push_back(cc); 324 } 325 } 326 } 327 } 328 return ans; 329 } 330 331 // TODO(sggonzalez): Adapt this function for the different color spaces or 332 // remove it if the color space with the best performance does not need it 333 void ComputeDotLosses(GaussianEllipse* ellipse, const ConnectedComponent& cc, 334 const Rect& rect, const Image3F& img, 335 const Image3F& background) { 336 const int rectBounds = 2; 337 const double kIntensityR = 0.0; // 0.015; 338 const double kSigmaR = 0.0; // 0.01; 339 const double kZeroEpsilon = 0.1; // Tolerance to consider a value negative 340 double ct = cos(ellipse->angle); 341 double st = sin(ellipse->angle); 342 const std::array<double, 3> channelGains{{1.0, 1.0, 1.0}}; 343 int N = 0; 344 ellipse->l1_loss = 0.0; 345 ellipse->l2_loss = 0.0; 346 ellipse->neg_pixels = 0; 347 ellipse->neg_value.fill(0.0); 348 double distMeanModeSq = (cc.mode.x - ellipse->x) * (cc.mode.x - ellipse->x) + 349 (cc.mode.y - ellipse->y) * (cc.mode.y - ellipse->y); 350 ellipse->custom_loss = 0.0; 351 for (int c = 0; c < 3; c++) { 352 for (int sy = -rectBounds; 353 sy < (static_cast<int>(cc.bounds.ysize()) + rectBounds); sy++) { 354 int y = sy + cc.bounds.y0(); 355 if (y < 0 || static_cast<size_t>(y) >= rect.ysize()) continue; 356 const float* JXL_RESTRICT row = rect.ConstPlaneRow(img, c, y); 357 // bgrow is only used if kOptimizeBackground is false. 358 // NOLINTNEXTLINE(clang-analyzer-deadcode.DeadStores) 359 const float* JXL_RESTRICT bgrow = rect.ConstPlaneRow(background, c, y); 360 for (int sx = -rectBounds; 361 sx < (static_cast<int>(cc.bounds.xsize()) + rectBounds); sx++) { 362 int x = sx + cc.bounds.x0(); 363 if (x < 0 || static_cast<size_t>(x) >= rect.xsize()) continue; 364 double target = row[x]; 365 double dotDelta = DotGaussianModel( 366 x - ellipse->x, y - ellipse->y, ct, st, ellipse->sigma_x, 367 ellipse->sigma_y, ellipse->intensity[c]); 368 if (dotDelta > target + kZeroEpsilon) { 369 ellipse->neg_pixels++; 370 ellipse->neg_value[c] += dotDelta - target; 371 } 372 double bkg = kOptimizeBackground ? ellipse->bgColor[c] : bgrow[x]; 373 double pred = bkg + dotDelta; 374 double diff = target - pred; 375 double l2 = channelGains[c] * diff * diff; 376 double l1 = channelGains[c] * std::fabs(diff); 377 ellipse->l2_loss += l2; 378 ellipse->l1_loss += l1; 379 double w = DotGaussianModel(x - cc.mode.x, y - cc.mode.y, 1.0, 0.0, 380 1.0 + ellipse->sigma_x, 381 1.0 + ellipse->sigma_y, 1.0); 382 ellipse->custom_loss += w * l2; 383 N++; 384 } 385 } 386 } 387 ellipse->l2_loss /= N; 388 ellipse->custom_loss /= N; 389 ellipse->custom_loss += 20.0 * distMeanModeSq + ellipse->neg_value[1]; 390 ellipse->l1_loss /= N; 391 double ridgeTerm = kSigmaR * ellipse->sigma_x + kSigmaR * ellipse->sigma_y; 392 for (int c = 0; c < 3; c++) { 393 ridgeTerm += kIntensityR * ellipse->intensity[c] * ellipse->intensity[c]; 394 } 395 ellipse->ridge_loss = ellipse->l2_loss + ridgeTerm; 396 } 397 398 GaussianEllipse FitGaussianFast(const ConnectedComponent& cc, const Rect& rect, 399 const Image3F& img, const Image3F& background) { 400 constexpr bool leastSqIntensity = true; 401 constexpr double kEpsilon = 1e-6; 402 GaussianEllipse ans; 403 constexpr int kRectBounds = (kEllipseWindowSize >> 1); 404 405 // Compute the 1st and 2nd moments of the CC 406 double sum = 0.0; 407 int N = 0; 408 std::array<double, 3> m1{{0.0, 0.0, 0.0}}; 409 std::array<double, 3> m2{{0.0, 0.0, 0.0}}; 410 std::array<double, 3> color{{0.0, 0.0, 0.0}}; 411 std::array<double, 3> bgColor{{0.0, 0.0, 0.0}}; 412 413 JXL_DEBUG(JXL_DEBUG_DOT_DETECT, 414 "%" PRIuS " %" PRIuS " %" PRIuS " %" PRIuS "\n", cc.bounds.x0(), 415 cc.bounds.y0(), cc.bounds.xsize(), cc.bounds.ysize()); 416 for (int c = 0; c < 3; c++) { 417 color[c] = rect.ConstPlaneRow(img, c, cc.mode.y)[cc.mode.x] - 418 rect.ConstPlaneRow(background, c, cc.mode.y)[cc.mode.x]; 419 } 420 double sign = (color[1] > 0) ? 1 : -1; 421 for (int sy = -kRectBounds; sy <= kRectBounds; sy++) { 422 int y = sy + cc.mode.y; 423 if (y < 0 || static_cast<size_t>(y) >= rect.ysize()) continue; 424 const float* JXL_RESTRICT row = rect.ConstPlaneRow(img, 1, y); 425 const float* JXL_RESTRICT bgrow = rect.ConstPlaneRow(background, 1, y); 426 for (int sx = -kRectBounds; sx <= kRectBounds; sx++) { 427 int x = sx + cc.mode.x; 428 if (x < 0 || static_cast<size_t>(x) >= rect.xsize()) continue; 429 double w = std::max(kEpsilon, sign * (row[x] - bgrow[x])); 430 sum += w; 431 432 m1[0] += w * x; 433 m1[1] += w * y; 434 m2[0] += w * x * x; 435 m2[1] += w * x * y; 436 m2[2] += w * y * y; 437 for (int c = 0; c < 3; c++) { 438 bgColor[c] += rect.ConstPlaneRow(background, c, y)[x]; 439 } 440 N++; 441 } 442 } 443 JXL_CHECK(N > 0); 444 445 for (int i = 0; i < 3; i++) { 446 m1[i] /= sum; 447 m2[i] /= sum; 448 bgColor[i] /= N; 449 } 450 451 // Some magic constants 452 constexpr double kSigmaMult = 1.0; 453 constexpr std::array<double, 3> kScaleMult{{1.1, 1.1, 1.1}}; 454 455 // Now set the parameters of the Gaussian 456 ans.x = m1[0]; 457 ans.y = m1[1]; 458 for (int j = 0; j < 3; j++) { 459 ans.intensity[j] = kScaleMult[j] * color[j]; 460 } 461 462 Matrix2x2 Sigma; 463 Vector2 d; 464 Matrix2x2 U; 465 Sigma[0][0] = m2[0] - m1[0] * m1[0]; 466 Sigma[1][1] = m2[2] - m1[1] * m1[1]; 467 Sigma[0][1] = Sigma[1][0] = m2[1] - m1[0] * m1[1]; 468 ConvertToDiagonal(Sigma, d, U); 469 Vector2& u = U[1]; 470 int p1 = 0; 471 int p2 = 1; 472 if (d[0] < d[1]) std::swap(p1, p2); 473 ans.sigma_x = kSigmaMult * d[p1]; 474 ans.sigma_y = kSigmaMult * d[p2]; 475 ans.angle = std::atan2(u[p1], u[p2]); 476 ans.l2_loss = 0.0; 477 ans.bgColor = bgColor; 478 if (leastSqIntensity) { 479 GaussianEllipse* ellipse = &ans; 480 double ct = cos(ans.angle); 481 double st = sin(ans.angle); 482 // Estimate intensity with least squares (fixed background) 483 for (int c = 0; c < 3; c++) { 484 double gg = 0.0; 485 double gd = 0.0; 486 int yc = static_cast<int>(cc.mode.y); 487 int xc = static_cast<int>(cc.mode.x); 488 for (int y = yc - kRectBounds; y <= yc + kRectBounds; y++) { 489 if (y < 0 || static_cast<size_t>(y) >= rect.ysize()) continue; 490 const float* JXL_RESTRICT row = rect.ConstPlaneRow(img, c, y); 491 const float* JXL_RESTRICT bgrow = rect.ConstPlaneRow(background, c, y); 492 for (int x = xc - kRectBounds; x <= xc + kRectBounds; x++) { 493 if (x < 0 || static_cast<size_t>(x) >= rect.xsize()) continue; 494 double target = row[x] - bgrow[x]; 495 double gaussian = 496 DotGaussianModel(x - ellipse->x, y - ellipse->y, ct, st, 497 ellipse->sigma_x, ellipse->sigma_y, 1.0); 498 gg += gaussian * gaussian; 499 gd += gaussian * target; 500 } 501 } 502 ans.intensity[c] = gd / (gg + 1e-6); // Regularized least squares 503 } 504 } 505 ComputeDotLosses(&ans, cc, rect, img, background); 506 return ans; 507 } 508 509 GaussianEllipse FitGaussian(const ConnectedComponent& cc, const Rect& rect, 510 const Image3F& img, const Image3F& background) { 511 auto ellipse = FitGaussianFast(cc, rect, img, background); 512 if (ellipse.sigma_x < ellipse.sigma_y) { 513 std::swap(ellipse.sigma_x, ellipse.sigma_y); 514 ellipse.angle += kPi / 2.0; 515 } 516 ellipse.angle -= kPi * std::floor(ellipse.angle / kPi); 517 if (fabs(ellipse.angle - kPi) < 1e-6 || fabs(ellipse.angle) < 1e-6) { 518 ellipse.angle = 0.0; 519 } 520 JXL_CHECK(ellipse.angle >= 0 && ellipse.angle <= kPi && 521 ellipse.sigma_x >= ellipse.sigma_y); 522 JXL_DEBUG(JXL_DEBUG_DOT_DETECT, 523 "Ellipse mu=(%lf,%lf) sigma=(%lf,%lf) angle=%lf " 524 "intensity=(%lf,%lf,%lf) bg=(%lf,%lf,%lf) l2_loss=%lf " 525 "custom_loss=%lf, neg_pix=%" PRIuS ", neg_v=(%lf,%lf,%lf)\n", 526 ellipse.x, ellipse.y, ellipse.sigma_x, ellipse.sigma_y, 527 ellipse.angle, ellipse.intensity[0], ellipse.intensity[1], 528 ellipse.intensity[2], ellipse.bgColor[0], ellipse.bgColor[1], 529 ellipse.bgColor[2], ellipse.l2_loss, ellipse.custom_loss, 530 ellipse.neg_pixels, ellipse.neg_value[0], ellipse.neg_value[1], 531 ellipse.neg_value[2]); 532 return ellipse; 533 } 534 535 } // namespace 536 537 StatusOr<std::vector<PatchInfo>> DetectGaussianEllipses( 538 const Image3F& opsin, const Rect& rect, const GaussianDetectParams& params, 539 const EllipseQuantParams& qParams, ThreadPool* pool) { 540 std::vector<PatchInfo> dots; 541 JXL_ASSIGN_OR_RETURN(Image3F smooth, 542 Image3F::Create(opsin.xsize(), opsin.ysize())); 543 JXL_ASSIGN_OR_RETURN(ImageF energy, ComputeEnergyImage(opsin, &smooth, pool)); 544 JXL_ASSIGN_OR_RETURN(std::vector<ConnectedComponent> components, 545 FindCC(energy, rect, params.t_low, params.t_high, 546 params.maxWinSize, params.minScore)); 547 size_t numCC = 548 std::min(params.maxCC, (components.size() * params.percCC) / 100); 549 if (components.size() > numCC) { 550 std::sort( 551 components.begin(), components.end(), 552 [](const ConnectedComponent& a, const ConnectedComponent& b) -> bool { 553 return a.score > b.score; 554 }); 555 components.erase(components.begin() + numCC, components.end()); 556 } 557 for (const auto& cc : components) { 558 GaussianEllipse ellipse = FitGaussian(cc, rect, opsin, smooth); 559 if (ellipse.x < 0.0 || 560 std::ceil(ellipse.x) >= static_cast<double>(rect.xsize()) || 561 ellipse.y < 0.0 || 562 std::ceil(ellipse.y) >= static_cast<double>(rect.ysize())) { 563 continue; 564 } 565 if (ellipse.neg_pixels > params.maxNegPixels) continue; 566 double intensity = 0.21 * ellipse.intensity[0] + 567 0.72 * ellipse.intensity[1] + 568 0.07 * ellipse.intensity[2]; 569 double intensitySq = intensity * intensity; 570 // for (int c = 0; c < 3; c++) { 571 // intensitySq += ellipse.intensity[c] * ellipse.intensity[c]; 572 //} 573 double sqDistMeanMode = (ellipse.x - cc.mode.x) * (ellipse.x - cc.mode.x) + 574 (ellipse.y - cc.mode.y) * (ellipse.y - cc.mode.y); 575 if (ellipse.l2_loss < params.maxL2Loss && 576 ellipse.custom_loss < params.maxCustomLoss && 577 intensitySq > (params.minIntensity * params.minIntensity) && 578 sqDistMeanMode < params.maxDistMeanMode * params.maxDistMeanMode) { 579 size_t x0 = cc.bounds.x0(); 580 size_t y0 = cc.bounds.y0(); 581 dots.emplace_back(); 582 dots.back().second.emplace_back(x0, y0); 583 QuantizedPatch& patch = dots.back().first; 584 patch.xsize = cc.bounds.xsize(); 585 patch.ysize = cc.bounds.ysize(); 586 for (size_t y = 0; y < patch.ysize; y++) { 587 for (size_t x = 0; x < patch.xsize; x++) { 588 for (size_t c = 0; c < 3; c++) { 589 patch.fpixels[c][y * patch.xsize + x] = 590 rect.ConstPlaneRow(opsin, c, y0 + y)[x0 + x] - 591 rect.ConstPlaneRow(smooth, c, y0 + y)[x0 + x]; 592 } 593 } 594 } 595 } 596 } 597 return dots; 598 } 599 600 } // namespace jxl 601 #endif // HWY_ONCE