|
9 | 9 | import os |
10 | 10 | import unittest |
11 | 11 | from common import run_tests |
12 | | -from tempfile import mktemp |
| 12 | +from numpy.linalg import eigh |
13 | 13 |
|
| 14 | +import math |
14 | 15 | import numpy as np |
15 | 16 | import numpy.testing as nptest |
16 | 17 | import pandas as pd |
| 18 | +import tempfile |
17 | 19 |
|
18 | 20 | import ants |
19 | 21 |
|
@@ -524,5 +526,72 @@ def test_label_image_registration(self): |
524 | 526 | moving_intensity_images=mi) |
525 | 527 |
|
526 | 528 |
|
| 529 | + def test_tensor_reorient(self): |
| 530 | + fa = 0.75 |
| 531 | + trace = 2.1E-3 |
| 532 | + l_1 = trace * (1 + (2 * fa) / np.sqrt(3 - 2 * fa**2)) |
| 533 | + l_2 = (trace - l_1) / 2 |
| 534 | + |
| 535 | + tensor = np.zeros((3, 3)) |
| 536 | + tensor[0, 0] = l_1 |
| 537 | + tensor[1, 1] = l_2 |
| 538 | + tensor[2, 2] = l_2 |
| 539 | + |
| 540 | + img_shape = (10,10,10) |
| 541 | + |
| 542 | + # Fill image with this tensor (SymmetricSecondRankTensor has 6 unique values) |
| 543 | + image_data = np.zeros(img_shape + (6,), dtype=np.float32) |
| 544 | + image_data[..., 0] = tensor[0, 0] # xx |
| 545 | + image_data[..., 1] = tensor[0, 1] # xy |
| 546 | + image_data[..., 2] = tensor[0, 2] # xz |
| 547 | + image_data[..., 3] = tensor[1, 1] # yy |
| 548 | + image_data[..., 4] = tensor[1, 2] # yz |
| 549 | + image_data[..., 5] = tensor[2, 2] # zz |
| 550 | + |
| 551 | + tensor_image = ants.from_numpy(image_data, has_components=True) |
| 552 | + tensor_image.set_direction(np.eye(3)) |
| 553 | + tensor_image.set_spacing((1.0, 1.0, 1.0)) |
| 554 | + tensor_image.set_origin((0.0, 0.0, 0.0)) |
| 555 | + |
| 556 | + ref_image = ants.from_numpy(image_data[..., 0], has_components=False) |
| 557 | + ants.copy_image_info(tensor_image, ref_image) |
| 558 | + |
| 559 | + theta = math.radians(20) |
| 560 | + tx_params = np.array([math.cos(theta), -math.sin(theta), 0, math.sin(theta), math.cos(theta), 0, 0, 0, 1, 0, 0, 0]) |
| 561 | + |
| 562 | + affine_tx = ants.create_ants_transform(transform_type="AffineTransform", dimension=3, parameters=tx_params) |
| 563 | + |
| 564 | + fd, tx_fn = tempfile.mkstemp(suffix=".mat") |
| 565 | + os.close(fd) |
| 566 | + |
| 567 | + ants.write_transform(affine_tx, tx_fn) |
| 568 | + |
| 569 | + reoriented_tensor = ants.apply_transforms(ref_image, tensor_image, tx_fn, imagetype=2, verbose=True).numpy() |
| 570 | + |
| 571 | + os.remove(tx_fn) |
| 572 | + |
| 573 | + mid = reoriented_tensor.shape[0] // 2 |
| 574 | + |
| 575 | + dtUpper = reoriented_tensor[mid, mid, mid] |
| 576 | + |
| 577 | + # Convert upper-triangular to full tensor |
| 578 | + T = np.array([ |
| 579 | + [dtUpper[0], dtUpper[1], dtUpper[2]], |
| 580 | + [dtUpper[1], dtUpper[3], dtUpper[4]], |
| 581 | + [dtUpper[2], dtUpper[4], dtUpper[5]], |
| 582 | + ]) |
| 583 | + |
| 584 | + # Check tensor has been reoriented by theta |
| 585 | + # Compute eigenvectors and eigenvalues |
| 586 | + eigvals, eigvecs = eigh(T) |
| 587 | + principal_ev = eigvecs[:, np.argmax(eigvals)] |
| 588 | + |
| 589 | + # The expected reorientation is the opposite to the forward transform |
| 590 | + expected_ev = np.array([math.cos(-theta), math.sin(-theta), 0]) |
| 591 | + |
| 592 | + dot = np.abs(np.dot(principal_ev, expected_ev)) |
| 593 | + assert math.isclose(dot, 1.0, abs_tol=1e-3), f"Principal eigenvector {principal_ev}, expected {expected_ev}" |
| 594 | + |
| 595 | + |
527 | 596 | if __name__ == "__main__": |
528 | 597 | run_tests() |
0 commit comments