@@ -628,6 +628,208 @@ class CV_EXPORTS_W BinaryDescriptor : public Algorithm
628
628
629
629
cv::Mat_<float > tempVecLineFit; // the vector used in line fit function;
630
630
631
+ /* * Compare doubles by relative error.
632
+ The resulting rounding error after floating point computations
633
+ depend on the specific operations done. The same number computed by
634
+ different algorithms could present different rounding errors. For a
635
+ useful comparison, an estimation of the relative rounding error
636
+ should be considered and compared to a factor times EPS. The factor
637
+ should be related to the accumulated rounding error in the chain of
638
+ computation. Here, as a simplification, a fixed factor is used.
639
+ */
640
+ static int double_equal ( double a, double b )
641
+ {
642
+ double abs_diff, aa, bb, abs_max;
643
+ /* trivial case */
644
+ if ( a == b )
645
+ return true ;
646
+ abs_diff = fabs ( a - b );
647
+ aa = fabs ( a );
648
+ bb = fabs ( b );
649
+ abs_max = aa > bb ? aa : bb;
650
+
651
+ /* DBL_MIN is the smallest normalized number, thus, the smallest
652
+ number whose relative error is bounded by DBL_EPSILON. For
653
+ smaller numbers, the same quantization steps as for DBL_MIN
654
+ are used. Then, for smaller numbers, a meaningful "relative"
655
+ error should be computed by dividing the difference by DBL_MIN. */
656
+ if ( abs_max < DBL_MIN )
657
+ abs_max = DBL_MIN;
658
+
659
+ /* equal if relative error <= factor x eps */
660
+ return ( abs_diff / abs_max ) <= ( RELATIVE_ERROR_FACTOR * DBL_EPSILON );
661
+ }
662
+
663
+ /* * Computes the natural logarithm of the absolute value of
664
+ the gamma function of x using the Lanczos approximation.
665
+ See http://www.rskey.org/gamma.htm
666
+ The formula used is
667
+ @f[
668
+ \Gamma(x) = \frac{ \sum_{n=0}^{N} q_n x^n }{ \Pi_{n=0}^{N} (x+n) }
669
+ (x+5.5)^{x+0.5} e^{-(x+5.5)}
670
+ @f]
671
+ so
672
+ @f[
673
+ \log\Gamma(x) = \log\left( \sum_{n=0}^{N} q_n x^n \right)
674
+ + (x+0.5) \log(x+5.5) - (x+5.5) - \sum_{n=0}^{N} \log(x+n)
675
+ @f]
676
+ and
677
+ q0 = 75122.6331530,
678
+ q1 = 80916.6278952,
679
+ q2 = 36308.2951477,
680
+ q3 = 8687.24529705,
681
+ q4 = 1168.92649479,
682
+ q5 = 83.8676043424,
683
+ q6 = 2.50662827511.
684
+ */
685
+ static double log_gamma_lanczos ( double x )
686
+ {
687
+ static double q[7 ] =
688
+ { 75122.6331530 , 80916.6278952 , 36308.2951477 , 8687.24529705 , 1168.92649479 , 83.8676043424 , 2.50662827511 };
689
+ double a = ( x + 0.5 ) * log ( x + 5.5 ) - ( x + 5.5 );
690
+ double b = 0.0 ;
691
+ int n;
692
+ for ( n = 0 ; n < 7 ; n++ )
693
+ {
694
+ a -= log ( x + (double ) n );
695
+ b += q[n] * pow ( x, (double ) n );
696
+ }
697
+ return a + log ( b );
698
+ }
699
+
700
+ /* * Computes the natural logarithm of the absolute value of
701
+ the gamma function of x using Windschitl method.
702
+ See http://www.rskey.org/gamma.htm
703
+ The formula used is
704
+ @f[
705
+ \Gamma(x) = \sqrt{\frac{2\pi}{x}} \left( \frac{x}{e}
706
+ \sqrt{ x\sinh(1/x) + \frac{1}{810x^6} } \right)^x
707
+ @f]
708
+ so
709
+ @f[
710
+ \log\Gamma(x) = 0.5\log(2\pi) + (x-0.5)\log(x) - x
711
+ + 0.5x\log\left( x\sinh(1/x) + \frac{1}{810x^6} \right).
712
+ @f]
713
+ This formula is a good approximation when x > 15.
714
+ */
715
+ static double log_gamma_windschitl ( double x )
716
+ {
717
+ return 0.918938533204673 + ( x - 0.5 ) * log ( x ) - x + 0.5 * x * log ( x * sinh ( 1 / x ) + 1 / ( 810.0 * pow ( x, 6.0 ) ) );
718
+ }
719
+
720
+ /* * Computes -log10(NFA).
721
+ NFA stands for Number of False Alarms:
722
+ @f[
723
+ \mathrm{NFA} = NT \cdot B(n,k,p)
724
+ @f]
725
+ - NT - number of tests
726
+ - B(n,k,p) - tail of binomial distribution with parameters n,k and p:
727
+ @f[
728
+ B(n,k,p) = \sum_{j=k}^n
729
+ \left(\begin{array}{c}n\\j\end{array}\right)
730
+ p^{j} (1-p)^{n-j}
731
+ @f]
732
+ The value -log10(NFA) is equivalent but more intuitive than NFA:
733
+ - -1 corresponds to 10 mean false alarms
734
+ - 0 corresponds to 1 mean false alarm
735
+ - 1 corresponds to 0.1 mean false alarms
736
+ - 2 corresponds to 0.01 mean false alarms
737
+ - ...
738
+ Used this way, the bigger the value, better the detection,
739
+ and a logarithmic scale is used.
740
+ @param n,k,p binomial parameters.
741
+ @param logNT logarithm of Number of Tests
742
+ The computation is based in the gamma function by the following
743
+ relation:
744
+ @f[
745
+ \left(\begin{array}{c}n\\k\end{array}\right)
746
+ = \frac{ \Gamma(n+1) }{ \Gamma(k+1) \cdot \Gamma(n-k+1) }.
747
+ @f]
748
+ We use efficient algorithms to compute the logarithm of
749
+ the gamma function.
750
+ To make the computation faster, not all the sum is computed, part
751
+ of the terms are neglected based on a bound to the error obtained
752
+ (an error of 10% in the result is accepted).
753
+ */
754
+ static double nfa ( int n, int k, double p, double logNT )
755
+ {
756
+ double tolerance = 0.1 ; /* an error of 10% in the result is accepted */
757
+ double log1term, term, bin_term, mult_term, bin_tail, err, p_term;
758
+ int i;
759
+
760
+ /* check parameters */
761
+ if ( n < 0 || k < 0 || k > n || p <= 0.0 || p >= 1.0 )
762
+ CV_Error (Error::StsBadArg, " nfa: wrong n, k or p values.\n " );
763
+ /* trivial cases */
764
+ if ( n == 0 || k == 0 )
765
+ return -logNT;
766
+ if ( n == k )
767
+ return -logNT - (double ) n * log10 ( p );
768
+
769
+ /* probability term */
770
+ p_term = p / ( 1.0 - p );
771
+
772
+ /* compute the first term of the series */
773
+ /*
774
+ binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
775
+ where bincoef(n,i) are the binomial coefficients.
776
+ But
777
+ bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
778
+ We use this to compute the first term. Actually the log of it.
779
+ */
780
+ log1term = log_gamma ( (double ) n + 1.0 )- log_gamma ( (double ) k + 1.0 )- log_gamma ( (double ) ( n - k ) + 1.0 )
781
+ + (double ) k * log ( p )
782
+ + (double ) ( n - k ) * log ( 1.0 - p );
783
+ term = exp ( log1term );
784
+
785
+ /* in some cases no more computations are needed */
786
+ if ( double_equal ( term, 0.0 ) )
787
+ { /* the first term is almost zero */
788
+ if ( (double ) k > (double ) n * p ) /* at begin or end of the tail? */
789
+ return -log1term / MLN10 - logNT; /* end: use just the first term */
790
+ else
791
+ return -logNT; /* begin: the tail is roughly 1 */
792
+ }
793
+
794
+ /* compute more terms if needed */
795
+ bin_tail = term;
796
+ for ( i = k + 1 ; i <= n; i++ )
797
+ {
798
+ /* As
799
+ term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
800
+ and
801
+ bincoef(n,i)/bincoef(n,i-1) = n-i+1 / i,
802
+ then,
803
+ term_i / term_i-1 = (n-i+1)/i * p/(1-p)
804
+ and
805
+ term_i = term_i-1 * (n-i+1)/i * p/(1-p).
806
+ p/(1-p) is computed only once and stored in 'p_term'.
807
+ */
808
+ bin_term = (double ) ( n - i + 1 ) / (double ) i;
809
+ mult_term = bin_term * p_term;
810
+ term *= mult_term;
811
+ bin_tail += term;
812
+ if ( bin_term < 1.0 )
813
+ {
814
+ /* When bin_term<1 then mult_term_j<mult_term_i for j>i.
815
+ Then, the error on the binomial tail when truncated at
816
+ the i term can be bounded by a geometric series of form
817
+ term_i * sum mult_term_i^j. */
818
+ err = term * ( ( 1.0 - pow ( mult_term, (double ) ( n - i + 1 ) ) ) / ( 1.0 - mult_term ) - 1.0 );
819
+ /* One wants an error at most of tolerance*final_result, or:
820
+ tolerance * abs(-log10(bin_tail)-logNT).
821
+ Now, the error that can be accepted on bin_tail is
822
+ given by tolerance*final_result divided by the derivative
823
+ of -log10(x) when x=bin_tail. that is:
824
+ tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
825
+ Finally, we truncate the tail if the error is less than:
826
+ tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */
827
+ if ( err < tolerance * fabs ( -log10 ( bin_tail ) - logNT ) * bin_tail )
828
+ break ;
829
+ }
830
+ }
831
+ return -log10 ( bin_tail ) - logNT;
832
+ }
631
833
};
632
834
633
835
// Specifies a vector of lines.
0 commit comments