| 901 | #define sign(x) (((x) < 0) ? (-1) : (1)) |
| 902 | |
| 903 | void vtkPKdTree::_select(int L, int R, int K, int dim) |
| 904 | { |
| 905 | int N, I, J, S, SD, LL, RR; |
| 906 | float Z; |
| 907 | |
| 908 | while (R > L) |
| 909 | { |
| 910 | if (R - L > 600) |
| 911 | { |
| 912 | // "Recurse on a sample of size S to get an estimate for the |
| 913 | // (K-L+1)-th smallest element into X[K], biased slightly so |
| 914 | // that the (K-L+1)-th element is expected to lie in the |
| 915 | // smaller set after partitioning" |
| 916 | |
| 917 | N = R - L + 1; |
| 918 | I = K - L + 1; |
| 919 | Z = static_cast<float>(log(float(N))); |
| 920 | S = static_cast<int>(.5 * exp(2 * Z / 3)); |
| 921 | SD = static_cast<int>(.5 * sqrt(Z * S * ((float)(N - S) / N)) * sign(I - N / 2)); |
| 922 | LL = vtkMath::Max(L, K - static_cast<int>((I * ((float)S / N))) + SD); |
| 923 | RR = vtkMath::Min(R, K + static_cast<int>((N - I) * ((float)S / N)) + SD); |
| 924 | this->_select(LL, RR, K, dim); |
| 925 | } |
| 926 | |
| 927 | int p1 = this->WhoHas(L); |
| 928 | int p2 = this->WhoHas(R); |
| 929 | |
| 930 | // "now adjust L,R so they surround the subset containing |
| 931 | // the (K-L+1)-th smallest element" |
| 932 | |
| 933 | // Due to very severe worst case behavior when the |
| 934 | // value at K (call it "T") is repeated many times in the array, we |
| 935 | // rearrange the array into three intervals: the leftmost being values |
| 936 | // less than T, the center being values equal to T, and the rightmost |
| 937 | // being values greater than T. Two integers are returned. This first |
| 938 | // is the global index of the start of the second interval. The second |
| 939 | // is the global index of the start of the third interval. (If there |
| 940 | // are no values greater than "T", the second integer will be R+1.) |
| 941 | // |
| 942 | // The original Floyd&Rivest arranged the array into two intervals, |
| 943 | // one less than "T", one greater than (or equal to) "T". |
| 944 | |
| 945 | int* idx = this->PartitionSubArray(L, R, K, dim, p1, p2); |
| 946 | |
| 947 | I = idx[0]; |
| 948 | J = idx[1]; |
| 949 | |
| 950 | if (K >= J) |
| 951 | { |
| 952 | L = J; |
| 953 | } |
| 954 | else if (K >= I) |
| 955 | { |
| 956 | L = R; // partitioning is done, K is in the interval of T's |
| 957 | } |
| 958 | else |
| 959 | { |
| 960 | R = I - 1; |