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