<jupyter_start><jupyter_text>Test of basic operations  rotate, keep simple surfaces, keep only bottom parts<jupyter_code>import sys
import os
import glob
import numpy as np
import math
from scipy import ndimage
from soma import aims
import anatomist.notebook as ana
from scipy.spatial.transform import Rotation as R

print(sys.version)<jupyter_output>3.6.9 (default, Jan 26 2021, 15:33:00) 
[GCC 8.4.0]<jupyter_text>Rotation with numpy / scipy<jupyter_code># This directory contains STS branch crops
crop_dir = '/neurospin/dico/data/deep_folding/data/crops/STS_branches/nearest/original/Lcrops'
##### Takes the first nii.gz file from crop_dir
nii_file = glob.glob(os.path.join(crop_dir, "*.nii.gz"))[0]
print(nii_file)
cropped_skeleton = aims.read(nii_file)<jupyter_output><empty_output><jupyter_text>The scope below is to visualize the file with anatomist:<jupyter_code>a = ana.Anatomist()
mri = a.loadObject(nii_file)
w1 = a.createWindow('Axial')
w2 = a.createWindow('Coronal')
mri.addInWindows([w1, w2])
w1.moveLinkedCursor([35., 35., 35.])<jupyter_output><empty_output><jupyter_text>We will now transform the nii.gz file into a numpy array and visualize it with Anatomist:<jupyter_code>arr_skel = np.asarray(cropped_skeleton)
np.unique(arr_skel)
block = a.createWindowsBlock(2)
w1 = a.createWindow('Axial', block=block)
w2 = a.createWindow('Axial', block=block)
img = aims.Volume(arr_skel)
img.header()['voxel_size'] = [1, 1, 1]
a_img = a.toAObject(img)
a_img.releaseAppRef()
a_img.setName("toto")
a_img.setChanged()
a_img.notifyObservers()
w1.addObjects(mri)
w2.addObjects([a_img])
w1.moveLinkedCursor([35., 35., 35.])
arr_skel.shape
arr_90 = np.rot90(arr_skel, axes=(0,1))
arr_90.shape
np.unique(arr_90)
block = a.createWindowsBlock(2)
w1 = a.createWindow('Axial', block=block)
w2 = a.createWindow('Axial', block=block)
img_90 = aims.Volume(arr_90)
img_90.header()['voxel_size'] = [1, 1, 1]
a_img_90 = a.toAObject(img_90)
a_img_90.releaseAppRef()
a_img_90.setName("toto")
a_img_90.setChanged()
a_img_90.notifyObservers()
w1.addObjects(mri)
w2.addObjects([a_img_90])
w1.moveLinkedCursor([35., 35., 35.])
def rotation_3d(a1, a2, a3):
    '''
    Return matrix 4x4 of rotation along 3 canonic axis based on 3 angles.

    - a1, a2, a3 : rotation angles.
    '''
    c1, s1 = np.cos(a1), np.sin(a1)
    c2, s2 = np.cos(a2), np.sin(a2)
    c3, s3 = np.cos(a3), np.sin(a3)
    m1 = np.matrix([
        [c1, -s1, 0, 0],
        [s1, c1,  0, 0],
        [0,   0,  1, 0],
        [0,   0,  0, 1]])
    m2 = np.matrix([
        [1,  0,   0, 0],
        [0, c2, -s2, 0],
        [0, s2,  c2, 0],
        [0,  0,   0, 1]])
    m3 = np.matrix([
        [c3, 0, -s3, 0],
        [0,  1,   0, 0],
        [s3, 0,  c3, 0],
        [0,  0,   0, 1]])
    return m1 * m2 * m3
rotation = rotation_3d(0., np.pi/8, 0.)
# rotation = rotation_3d(0, 0, 0)
translate = np.array([18., 8., 4., 1.])
translation = rotation_3d(0., 0., 0.)
translation[:, 3] = np.asmatrix(translate).T

print(rotation, "\n\n", translation)
tr_rot = aims.AffineTransformation3d(rotation)
tr_trans = aims.AffineTransformation3d(translation)
#tr = aims.AffineTransformation3d(r)
#tr.rotationaroundx(math.pi/16)
#tr.setTranslation((18, 8, 4))
rf = getattr(aims, 'ResamplerFactory_%s' % aims.voxelTypeCode(cropped_skeleton))()
resampler = rf.getResampler(0)

print(cropped_skeleton.header()['voxel_size'])
cropped_skeleton.shape
out_vol = aims.Volume(80, 80, 80, 1, dtype=aims.voxelTypeCode(cropped_skeleton))
padded_skeleton = aims.Volume(80, 80, 80, 1, dtype=aims.voxelTypeCode(cropped_skeleton))
out_vol.header()['voxel_size'] = [1., 1., 1.]
padded_skeleton.header()['voxel_size'] = [1., 1., 1.]
resampler.resample(cropped_skeleton, tr_trans, 0, padded_skeleton) 
resampler.resample(padded_skeleton, tr_rot, 0, out_vol) 
out_vol.shape
padded_skeleton.shape
out_arr = out_vol.arraydata()
np.unique(out_arr)
w1 = a.createWindow('Axial')
w2 = a.createWindow('Axial')
a_out = a.toAObject(out_vol)
a_out.addInWindows(w1)
mri.addInWindows(w2)
w1.moveLinkedCursor([35., 35., 35.])
type(out_vol)
cropped_skeleton.header()
out_vol.header()
padded_skeleton.header()
import skimage
im = skimage.data.astronaut()
np.unique(im)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15904761904761905, 'fraction_numerical': 0.054285714285714284, 'mean_word_length': 4.7785419532324624, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '324191', 'license': 'MIT', 'n_tokens_mistral': 1876, 'n_tokens_neox': 1760, 'n_words': 353, 'path': 'notebooks/test_basic_operations.ipynb', 'repo_name': 'neurospin-projects/2021_jchavas_dicoFolding', 'size': '4200'}