Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61388170
KendallTauRankCorrelation.java
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, May 6, 08:46
Size
10 KB
Mime Type
text/x-java
Expires
Wed, May 8, 08:46 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17505092
Attached To
R6269 BIOP Run Macro
KendallTauRankCorrelation.java
View Options
/*-
* #%L
* Fiji's plugin for colocalization analysis.
* %%
* Copyright (C) 2009 - 2017 Fiji developers.
* %%
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3 of the
* License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program. If not, see
* <http://www.gnu.org/licenses/gpl-3.0.html>.
* #L%
*/
package
sc.fiji.coloc.algorithms
;
import
ij.IJ
;
import
java.util.Arrays
;
import
net.imglib2.PairIterator
;
import
net.imglib2.RandomAccessible
;
import
net.imglib2.RandomAccessibleInterval
;
import
net.imglib2.TwinCursor
;
import
net.imglib2.type.logic.BitType
;
import
net.imglib2.type.numeric.RealType
;
import
net.imglib2.view.Views
;
import
sc.fiji.coloc.gadgets.DataContainer
;
import
sc.fiji.coloc.results.ResultHandler
;
/**
* This algorithm calculates Kendall's Tau-b rank correlation coefficient
* <p>
* According to <a
* href="http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient">this
* article</a>, Tau-b (appropriate if multiple pairs share the same first, or
* second, value), the rank correlation of a set of pairs <tt>(x_1, y_1), ...,
* (x_n, y_n)</tt>:
* </p>
*
* <pre>
* Tau_B = (n_c - n_d) / sqrt( (n_0 - n_1) (n_0 - n_2) )
* </pre>
*
* where
*
* <pre>
* n_0 = n (n - 1) / 2
* n_1 = sum_i t_i (t_i - 1) / 2
* n_2 = sum_j u_j (u_j - 1) / 2
* n_c = #{ i, j; i != j && (x_i - x_j) * (y_i - y_j) > 0 },
* i.e. the number of pairs of pairs agreeing on the order of x and y, respectively
* n_d = #{ i, j: i != j && (x_i - x_j) * (y_i - y_j) < 0 },
* i.e. the number of pairs of pairs where x and y are ordered opposite of each other
* t_i = number of tied values in the i-th group of ties for the first quantity
* u_j = number of tied values in the j-th group of ties for the second quantity
* </pre>
*
* @author Johannes Schindelin
* @param <T>
*/
public
class
KendallTauRankCorrelation
<
T
extends
RealType
<
T
>>
extends
Algorithm
<
T
>
{
public
KendallTauRankCorrelation
()
{
super
(
"Kendall's Tau-b Rank Correlation"
);
}
private
double
tau
;
@Override
public
void
execute
(
DataContainer
<
T
>
container
)
throws
MissingPreconditionException
{
RandomAccessible
<
T
>
img1
=
container
.
getSourceImage1
();
RandomAccessible
<
T
>
img2
=
container
.
getSourceImage2
();
RandomAccessibleInterval
<
BitType
>
mask
=
container
.
getMask
();
TwinCursor
<
T
>
cursor
=
new
TwinCursor
<
T
>(
img1
.
randomAccess
(),
img2
.
randomAccess
(),
Views
.
iterable
(
mask
).
localizingCursor
());
tau
=
calculateMergeSort
(
cursor
);
}
public
static
<
T
extends
RealType
<
T
>>
double
calculateNaive
(
final
PairIterator
<
T
>
iterator
)
{
if
(!
iterator
.
hasNext
())
{
return
Double
.
NaN
;
}
// See http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient
int
n
=
0
,
max1
=
0
,
max2
=
0
,
max
=
255
;
int
[][]
histogram
=
new
int
[
max
+
1
][
max
+
1
];
while
(
iterator
.
hasNext
())
{
iterator
.
fwd
();
T
type1
=
iterator
.
getFirst
();
T
type2
=
iterator
.
getSecond
();
double
ch1
=
type1
.
getRealDouble
();
double
ch2
=
type2
.
getRealDouble
();
if
(
ch1
<
0
||
ch2
<
0
||
ch1
>
max
||
ch2
>
max
)
{
IJ
.
log
(
"Error: The current Kendall Tau implementation is limited to 8-bit data"
);
return
Double
.
NaN
;
}
n
++;
int
ch1Int
=
(
int
)
Math
.
round
(
ch1
);
int
ch2Int
=
(
int
)
Math
.
round
(
ch2
);
histogram
[
ch1Int
][
ch2Int
]++;
if
(
max1
<
ch1Int
)
{
max1
=
ch1Int
;
}
if
(
max2
<
ch2Int
)
{
max2
=
ch2Int
;
}
}
long
n0
=
n
*
(
long
)(
n
-
1
)
/
2
,
n1
=
0
,
n2
=
0
,
nc
=
0
,
nd
=
0
;
for
(
int
i1
=
0
;
i1
<=
max1
;
i1
++)
{
IJ
.
log
(
""
+
i1
+
"/"
+
max1
);
int
ch1
=
0
;
for
(
int
i2
=
0
;
i2
<=
max2
;
i2
++)
{
ch1
+=
histogram
[
i1
][
i2
];
int
count
=
histogram
[
i1
][
i2
];
for
(
int
j1
=
0
;
j1
<
i1
;
j1
++)
{
for
(
int
j2
=
0
;
j2
<
i2
;
j2
++)
{
nc
+=
count
*
histogram
[
j1
][
j2
];
}
}
for
(
int
j1
=
0
;
j1
<
i1
;
j1
++)
{
for
(
int
j2
=
i2
+
1
;
j2
<=
max2
;
j2
++)
{
nd
+=
count
*
histogram
[
j1
][
j2
];
}
}
}
n1
+=
ch1
*
(
long
)(
ch1
-
1
)
/
2
;
}
for
(
int
i2
=
0
;
i2
<=
max2
;
i2
++)
{
int
ch2
=
0
;
for
(
int
i1
=
0
;
i1
<=
max1
;
i1
++)
{
ch2
+=
histogram
[
i1
][
i2
];
}
n2
+=
ch2
*
(
long
)(
ch2
-
1
)
/
2
;
}
return
(
nc
-
nd
)
/
Math
.
sqrt
((
n0
-
n1
)
*
(
double
)(
n0
-
n2
));
}
private
static
<
T
extends
RealType
<
T
>>
double
[][]
getPairs
(
final
PairIterator
<
T
>
iterator
)
{
// TODO: it is ridiculous that this has to be counted all the time (i.e. in most if not all measurements!).
// We only need an upper bound to begin with, so even the number of pixels in the first channel would be enough!
int
capacity
=
0
;
while
(
iterator
.
hasNext
())
{
iterator
.
fwd
();
capacity
++;
}
double
[]
values1
=
new
double
[
capacity
];
double
[]
values2
=
new
double
[
capacity
];
iterator
.
reset
();
int
count
=
0
;
while
(
iterator
.
hasNext
())
{
iterator
.
fwd
();
values1
[
count
]
=
iterator
.
getFirst
().
getRealDouble
();
values2
[
count
]
=
iterator
.
getSecond
().
getRealDouble
();
count
++;
}
if
(
count
<
capacity
)
{
values1
=
Arrays
.
copyOf
(
values1
,
count
);
values2
=
Arrays
.
copyOf
(
values2
,
count
);
}
return
new
double
[][]
{
values1
,
values2
};
}
/**
* Calculate Tau-b efficiently.
* <p>
* This implementation is based on <a
* href="http://en.wikipedia.org/wiki/Kendall_tau_rank_correlation_coefficient#Algorithms">this
* description of the merge sort based way to calculate Tau-b</a>. This is
* supposed to be the method described in:
* </p>
* <blockquote>Knight, W. (1966). "A Computer Method for Calculating
* Kendall's Tau with Ungrouped Data". Journal of the American Statistical
* Association 61 (314): 436–439. doi:10.2307/2282833.</blockquote>
* <p>
* but since that article is not available as Open Access, it is
* unnecessarily hard to verify.
* </p>
*
* @param iterator the iterator of the pairs
* @return Tau-b
*/
public
static
<
T
extends
RealType
<
T
>>
double
calculateMergeSort
(
final
PairIterator
<
T
>
iterator
)
{
final
double
[][]
pairs
=
getPairs
(
iterator
);
final
double
[]
x
=
pairs
[
0
];
final
double
[]
y
=
pairs
[
1
];
final
int
n
=
x
.
length
;
int
[]
index
=
new
int
[
n
];
for
(
int
i
=
0
;
i
<
n
;
i
++)
{
index
[
i
]
=
i
;
}
// First sort by x as primary key, y as secondary one.
// We use IntroSort here because it is fast and in-place.
IntArraySorter
.
sort
(
index
,
new
IntComparator
()
{
@Override
public
int
compare
(
int
a
,
int
b
)
{
double
xa
=
x
[
a
],
ya
=
y
[
a
];
double
xb
=
x
[
b
],
yb
=
y
[
b
];
int
result
=
Double
.
compare
(
xa
,
xb
);
return
result
!=
0
?
result
:
Double
.
compare
(
ya
,
yb
);
}
});
// The trick is to count the ties of x (n1) and the joint ties of x and y (n3) now, while
// index is sorted with regards to x.
long
n0
=
n
*
(
long
)(
n
-
1
)
/
2
;
long
n1
=
0
,
n3
=
0
;
for
(
int
i
=
1
;
i
<
n
;
i
++)
{
double
x0
=
x
[
index
[
i
-
1
]];
if
(
x
[
index
[
i
]]
!=
x0
)
{
continue
;
}
double
y0
=
y
[
index
[
i
-
1
]];
int
i1
=
i
;
do
{
double
y1
=
y
[
index
[
i1
++]];
if
(
y1
==
y0
)
{
int
i2
=
i1
;
while
(
i1
<
n
&&
x
[
index
[
i1
]]
==
x0
&&
y
[
index
[
i1
]]
==
y0
)
{
i1
++;
}
n3
+=
(
i1
-
i2
+
2
)
*
(
long
)(
i1
-
i2
+
1
)
/
2
;
}
y0
=
y1
;
}
while
(
i1
<
n
&&
x
[
index
[
i1
]]
==
x0
);
n1
+=
(
i1
-
i
+
1
)
*
(
long
)(
i1
-
i
)
/
2
;
i
=
i1
;
}
// Now, let's perform that merge sort that also counts S, the number of
// swaps a Bubble Sort would require (and which therefore is half the number
// by which we have to adjust n_0 - n_1 - n_2 + n_3 to obtain n_c - n_d)
final
MergeSort
mergeSort
=
new
MergeSort
(
index
,
new
IntComparator
()
{
@Override
public
int
compare
(
int
a
,
int
b
)
{
double
ya
=
y
[
a
];
double
yb
=
y
[
b
];
return
Double
.
compare
(
ya
,
yb
);
}
});
long
S
=
mergeSort
.
sort
();
index
=
mergeSort
.
getSorted
();
long
n2
=
0
;
for
(
int
i
=
1
;
i
<
n
;
i
++)
{
double
y0
=
y
[
index
[
i
-
1
]];
if
(
y
[
index
[
i
]]
!=
y0
)
{
continue
;
}
int
i1
=
i
+
1
;
while
(
i1
<
n
&&
y
[
index
[
i1
]]
==
y0
)
{
i1
++;
}
n2
+=
(
i1
-
i
+
1
)
*
(
long
)(
i1
-
i
)
/
2
;
i
=
i1
;
}
return
(
n0
-
n1
-
n2
+
n3
-
2
*
S
)
/
Math
.
sqrt
((
n0
-
n1
)
*
(
double
)(
n0
-
n2
));
}
private
final
static
class
MergeSort
{
private
int
[]
index
;
private
final
IntComparator
comparator
;
public
MergeSort
(
int
[]
index
,
IntComparator
comparator
)
{
this
.
index
=
index
;
this
.
comparator
=
comparator
;
}
public
int
[]
getSorted
()
{
return
index
;
}
/**
* Sorts the {@link #index} array.
* <p>
* This implements a non-recursive merge sort.
* </p>
* @param begin
* @param end
* @return the equivalent number of BubbleSort swaps
*/
public
long
sort
()
{
long
swaps
=
0
;
int
n
=
index
.
length
;
// There are merge sorts which perform in-place, but their runtime is worse than O(n log n)
int
[]
index2
=
new
int
[
n
];
for
(
int
step
=
1
;
step
<
n
;
step
<<=
1
)
{
int
begin
=
0
,
k
=
0
;
for
(;;)
{
int
begin2
=
begin
+
step
,
end
=
begin2
+
step
;
if
(
end
>=
n
)
{
if
(
begin2
>=
n
)
{
break
;
}
end
=
n
;
}
// calculate the equivalent number of BubbleSort swaps
// and perform merge, too
int
i
=
begin
,
j
=
begin2
;
while
(
i
<
begin2
&&
j
<
end
)
{
int
compare
=
comparator
.
compare
(
index
[
i
],
index
[
j
]);
if
(
compare
>
0
)
{
swaps
+=
(
begin2
-
i
);
index2
[
k
++]
=
index
[
j
++];
}
else
{
index2
[
k
++]
=
index
[
i
++];
}
}
if
(
i
<
begin2
)
{
do
{
index2
[
k
++]
=
index
[
i
++];
}
while
(
i
<
begin2
);
}
else
{
while
(
j
<
end
)
{
index2
[
k
++]
=
index
[
j
++];
}
}
begin
=
end
;
}
if
(
k
<
n
)
{
System
.
arraycopy
(
index
,
k
,
index2
,
k
,
n
-
k
);
}
int
[]
swapIndex
=
index2
;
index2
=
index
;
index
=
swapIndex
;
}
return
swaps
;
}
}
@Override
public
void
processResults
(
ResultHandler
<
T
>
handler
)
{
super
.
processResults
(
handler
);
handler
.
handleValue
(
"Kendall's Tau-b rank correlation value"
,
tau
,
4
);
}
}
Event Timeline
Log In to Comment