~nipy-developers/nipy/fff2 : revision 223.1.6

13

# Standard library imports

14

import os

15

import tempfile

16

import sys

16

17

18

# Standard scientific libraries imports (more specific imports are

18

19

# delayed, so that the part module can be used without them).

19

20

import numpy as np

20

21

import pylab as pl

22

from scipy import ndimage, stats

23

from nifti import NiftiImage

24

25

# Local imports

26

from fff2.utils.mask import compute_mask, _largest_cc

27

from fff2.utils.emp_null import ENN

28

import fff2.data

21

29

22

30

# The sform for MNI templates

23

31

mni_sform = np.array([[-1, 0, 0, 90],

39

47

# Functions for automatic choice of cuts coordinates

40

48

################################################################################

41

49

50

def largest_cc(mask):

51

""" Largest connect components extraction code that fails gracefully

52

when mask is too big for C extension code.

53

"""

54

try:

55

mask = _largest_cc(mask)

56

except TypeError:

57

pass

58

return mask.astype(np.bool)

59

42

60

def coord_transform(x, y, z, affine):

43

61

""" Convert the x, y, z coordinates from one image space to another

44

space, using the given affine, that maps from input to output

45

space.

46

47

Warning: The x, y and z have their Talairach meaning, not 3D

48

numy image meaning.

62

space.

63

64

Parameters

65

----------

66

x : number or ndarray

67

The x coordinates in the input space

68

y : number or ndarray

69

The y coordinates in the input space

70

z : number or ndarray

71

The z coordinates in the input space

72

affine : 2D 4x4 ndarray

73

affine that maps from input to output space.

74

75

Returns

76

-------

77

x : number or ndarray

78

The x coordinates in the output space

79

y : number or ndarray

80

The y coordinates in the output space

81

z : number or ndarray

82

The z coordinates in the output space

83

84

Warning: The x, y and z have their Talairach ordering, not 3D

85

numy image ordering.

49

86

"""

50

87

coords = np.c_[np.atleast_1d(x).flat,

51

88

np.atleast_1d(y).flat,

74

111

75

112

76

113

def find_activation(map, mask=None, pvalue=0.05, upper_only=False):

77

""" Use the fff2 FDR estimator to guess threshold for negative

78

and positive activation.

79

80

If mask is None, the mask is computed from the map.

81

If mask is False, the map is assumed to be already masked

82

(and thus flat).

83

If mask is a 3D boolean array, it is used to mask the map.

84

85

If upper_only is true, only a threshold for positive

86

activation is estimated.

114

""" Use the empirical normal null estimator to find threshold for

115

negative and positive activation.

116

117

Parameters

118

----------

119

map : 3D ndarray

120

The activation map, as a 3D image.

121

mask : 3D ndarray, boolean, optional

122

The brain mask. If None, the mask is computed from the map.

123

pvalue : float, optional

124

The pvalue of the false discovery rate used.

125

upper_only : boolean, optional

126

If true, only a threshold for positive activation is estimated.

127

128

Returns

129

-------

130

vmin : float, optional

131

The upper threshold for negative activation, not returned

132

if upper_only is True.

133

vmax : float

134

The lower threshod for positive activation.

87

135

"""

88

from fff2.utils.emp_null import ENN

89

from fff2.utils.mask import compute_mask_intra_array

90

136

91

137

if mask is None:

92

mask = compute_mask_intra_array(map)

93

if mask is not False:

94

map = map[mask]

138

mask = compute_mask(map)

139

map = map[mask]

95

140

vmax = ENN(map).threshold(alpha=pvalue)

96

141

if upper_only:

97

142

return vmax

101

146

102

147

def find_cut_coords(map, mask=None, activation_threshold=None):

103

148

""" Find the center of the largest activation connect component.

104

Returns the coordinates in voxels.

105

106

Mask is the brain mask.

107

108

If mask is None, the mask is computed from the map.

109

If mask is False, the map is assumed to be already masked

110

(and thus flat).

111

If mask is a 3D boolean array, it is used to mask the map.

149

150

Parameters

151

-----------

152

map : 3D ndarray

153

The activation map, as a 3D image.

154

mask : 3D ndarray, boolean, optional

155

The brain mask. If None, the mask is computed from the map.

156

activation_threshold : float, optional

157

The lower threshold to the positive activation. If None, the

158

activation threshold is computed using find_activation.

159

160

Returns

161

-------

162

x: float

163

the x coordinate in voxels.

164

y: float

165

the y coordinate in voxels.

166

z: float

167

the z coordinate in voxels.

112

168

"""

113

from fff2.utils.mask import _largest_cc

114

from scipy import ndimage, stats

115

169

#

116

170

my_map = map.copy()

117

171

if activation_threshold is None:

120

174

else:

121

175

mask = map>activation_threshold

122

176

if np.any(mask):

123

mask = _largest_cc(mask).astype(np.bool)

177

mask = largest_cc(mask)

124

178

my_map[np.logical_not(mask)] = 0

125

179

second_threshold = stats.scoreatpercentile(my_map[mask], 60)

126

180

if (my_map>second_threshold).sum() > 50:

127

my_map[np.logical_not(_largest_cc(my_map>second_threshold))] = 0

181

my_map[np.logical_not(largest_cc(my_map>second_threshold))] = 0

128

182

cut_coords = ndimage.center_of_mass(my_map)

129

183

return cut_coords

130

184

154

208

def get_anat(cls):

155

209

if cls.anat is not None:

156

210

return cls.anat, cls.anat_sform, cls.anat_max

157

# Do the import only here for loose coupling with pynifti and

158

# scipy

159

from scipy import ndimage

160

from nifti import NiftiImage

161

from fff2 import data

162

211

anat_im = NiftiImage(

163

212

os.path.join(os.path.dirname(

164

os.path.realpath(data.__file__)),

213

os.path.realpath(fff2.data.__file__)),

165

214

'MNI152_T1_1mm_brain.nii.gz'

166

215

))

167

216

anat = anat_im.data.T

177

226

def get_blurred(cls):

178

227

if cls.blurred is not None:

179

228

return cls.blurred

180

from scipy import ndimage

181

229

anat, _, _ = cls.get_anat()

182

230

cls.blurred = ndimage.gaussian_filter(

183

231

(ndimage.morphology.binary_fill_holes(

199

247

mask=None):

200

248

""" Plot three cuts of a given activation map (Frontal, Axial, and Lateral)

201

249

250

Parameters

251

----------

252

map : 3D ndarray

253

The activation map, as a 3D image.

254

sform : 4x4 ndarray

255

The affine matrix going from image voxel space to MNI space.

256

cut_coords: 3-tuple of floats

257

The MNI coordinates of the cut to perform, in MNI coordinates

258

and order.

259

anat : 3D ndarray, optional

260

The anatomical image to be used as a background. If None, the

261

MNI152 T1 1mm template is used.

262

anat_sform : 4x4 ndarray, optional

263

The affine matrix going from the anatomical image voxel space to

264

MNI space. This parameter is not used when the default

265

anatomical is used, but it is compulsory when using an

266

explicite anatomical image.

267

vmin : float, optional

268

The lower threshold of the positive activation. This

269

parameter is used to threshold the activation map.

270

figure_num : integer, optional

271

The number of the matplotlib figure used. If None is given, a

272

new figure is created.

273

axes : 4 tuple of float: (xmin, xmax, ymin, ymin), optional

274

The coordinates, in matplotlib figure space, of the axes

275

used to display the plot. If None, the complete figure is

276

used.

277

title : string, optional

278

The title dispayed on the figure.

279

mask : 3D ndarray, boolean, optional

280

The brain mask. If None, the mask is computed from the map.

281

282

Notes

283

-----

202

284

All the 3D arrays are in numpy convention: (x, y, z)

203

285

204

anat and anat_sform are respectively the 3D anatomical image used

205

for the background, and its sform. If None is given, a MNI

206

template is used.

207

208

286

Cut coordinates are in Talairach coordinates. Warning: Talairach

209

287

coordinates are (y, x, z), if (x, y, z) are in voxel-ordering

210

288

convention.

211

212

axes is an optional four-tuple [xmin, xmax, ymin, ymax] giving

213

the region, in figure units, of the figure used to plot.

214

289

"""

215

290

if anat is None:

216

291

anat, anat_sform, vmax_anat = _AnatCache.get_anat()

377

452

x, y, z = -6, -53, 9

378

453

x_map, y_map, z_map = coord_transform(x, y, z, mni_sform_inv)

379

454

map[x_map-30:x_map+30, y_map-3:y_map+3, z_map-10:z_map+10] = 1

380

map = map.T

381

455

map = np.ma.masked_less(map, 0.5)

382

456

plot_map_2d(map, mni_sform, cut_coords=(x, y, z),

383

457

figure_num=512)

388

462

""" Plot a 3D volume rendering view of the activation, with an

389

463

outline of the brain.

390

464

465

Parameters

466

----------

467

map : 3D ndarray

468

The activation map, as a 3D image.

469

sform : 4x4 ndarray

470

The affine matrix going from image voxel space to MNI space.

471

cut_coords: 3-tuple of floats, optional

472

The MNI coordinates of a 3D cursor to indicate a feature

473

or a cut, in MNI coordinates and order.

474

anat : 3D ndarray, optional

475

The anatomical image to be used as a background. If None, the

476

MNI152 T1 1mm template is used.

477

anat_sform : 4x4 ndarray, optional

478

The affine matrix going from the anatomical image voxel space to

479

MNI space. This parameter is not used when the default

480

anatomical is used, but it is compulsory when using an

481

explicite anatomical image.

482

vmin : float, optional

483

The lower threshold of the positive activation. This

484

parameter is used to threshold the activation map.

485

figure_num : integer, optional

486

The number of the Mayavi figure used. If None is given, a

487

new figure is created.

488

mask : 3D ndarray, boolean, optional

489

The brain mask. If None, the mask is computed from the map.

490

491

Notes

492

-----

391

493

All the 3D arrays are in numpy convention: (x, y, z)

392

494

393

anat and anat_sform are respectively the 3D anatomical image used

394

for the background, and its sform. If None is given, a MNI

395

template is used.

396

397

If cut_coords is given, it is Talairach coordinates of a

398

cross cursor position. Warning: Talairach coordinates are

399

(y, x, z), if (x, y, z) are in voxel-ordering convention.

400

401

The mask is used to estimate threshold for activation, if

402

necessary.

495

Cut coordinates are in Talairach coordinates. Warning: Talairach

496

coordinates are (y, x, z), if (x, y, z) are in voxel-ordering

497

convention.

498

499

If you are using a VTK version below 5.2, there is no way to

500

avoid opening a window during the rendering under Linux. This is

501

necessary to use the graphics card for the rendering. You must

502

maintain this window on top of others and on the screen.

503

504

If you are running on Windows, or using a recent version of VTK,

505

you can force offscreen mode using::

506

507

from enthought.mayavi import mlab

508

mlab.options.offscreen = True

403

509

"""

510

404

511

from enthought.mayavi import mlab

405

512

from enthought.mayavi.sources.api import ArraySource

406

513

if anat is None:

509

616

figure_num=512)

510

617

511

618

512

def plot_map(map, sform, cut_coords=None, anat=None, anat_sform=None,

619

def plot_map(map, sform, cut_coords, anat=None, anat_sform=None,

513

620

vmin=None, figure_num=None, title='', mask=None):

514

621

""" Plot a together a 3D volume rendering view of the activation, with an

515

622

outline of the brain, and 2D cuts. If Mayavi is not installed,

516

623

falls back to 2D views only.

517

624

625

Parameters

626

----------

627

map : 3D ndarray

628

The activation map, as a 3D image.

629

sform : 4x4 ndarray

630

The affine matrix going from image voxel space to MNI space.

631

cut_coords: 3-tuple of floats, optional

632

The MNI coordinates of the cut to perform, in MNI coordinates

633

and order. If None is given, the cut_coords are automaticaly

634

estimated.

635

anat : 3D ndarray, optional

636

The anatomical image to be used as a background. If None, the

637

MNI152 T1 1mm template is used.

638

anat_sform : 4x4 ndarray, optional

639

The affine matrix going from the anatomical image voxel space to

640

MNI space. This parameter is not used when the default

641

anatomical is used, but it is compulsory when using an

642

explicite anatomical image.

643

vmin : float, optional

644

The lower threshold of the positive activation. This

645

parameter is used to threshold the activation map.

646

figure_num : integer, optional

647

The number of the matplotlib and Mayavi figures used. If None is

648

given, a new figure is created.

649

title : string, optional

650

The title dispayed on the figure.

651

mask : 3D ndarray, boolean, optional

652

The brain mask. If None, the mask is computed from the map.

653

654

Notes

655

-----

518

656

All the 3D arrays are in numpy convention: (x, y, z)

519

657

520

anat and anat_sform are respectively the 3D anatomical image used

521

for the background, and its sform. If None is given, a MNI

522

template is used.

523

524

If cut_coords is given, it is Talairach coordinates of a

525

cross cursor position. Warning: Talairach coordinates are

526

(y, x, z), if (x, y, z) are in voxel-ordering convention.

658

Cut coordinates are in Talairach coordinates. Warning: Talairach

659

coordinates are (y, x, z), if (x, y, z) are in voxel-ordering

660

convention.

527

661

"""

528

662

try:

529

663

from enthought.mayavi import version

530

664

if not int(version.version[0]) > 2:

531

665

raise ImportError

532

666

except ImportError:

533

import sys

534

667

print >> sys.stderr, 'Mayavi > 3.x not installed, plotting only 2D'

535

668

return plot_map_2d(map, sform, cut_coords=cut_coords, anat=anat,

536

669

anat_sform=anat_sform, vmin=vmin,

562

695

x, y, z = -6, -53, 9

563

696

x_map, y_map, z_map = coord_transform(x, y, z, mni_sform_inv)

564

697

map[x_map-30:x_map+30, y_map-3:y_map+3, z_map-10:z_map+10] = 1

565

map = map.T

566

698

plot_map(map, mni_sform, cut_coords=(x, y, z), vmin=0.5,

567

699

figure_num=512)

568

700

570

702

def auto_plot_map(map, sform, vmin=None, cut_coords=None, do3d=False,

571

703

anat=None, anat_sform=None, title='',

572

704

figure_num=None, mask=None, auto_sign=True):

705

""" Automatic plotting of an activation map.

706

707

Plot a together a 3D volume rendering view of the activation, with an

708

outline of the brain, and 2D cuts. If Mayavi is not installed,

709

falls back to 2D views only.

710

711

Parameters

712

----------

713

map : 3D ndarray

714

The activation map, as a 3D image.

715

sform : 4x4 ndarray

716

The affine matrix going from image voxel space to MNI space.

717

vmin : float, optional

718

The lower threshold of the positive activation. This

719

parameter is used to threshold the activation map.

720

cut_coords: 3-tuple of floats, optional

721

The MNI coordinates of the cut to perform, in MNI coordinates

722

and order. If None is given, the cut_coords are automaticaly

723

estimated.

724

do3d : boolean, optional

725

If do3d is True, a 3D plot is created if Mayavi is installed.

726

anat : 3D ndarray, optional

727

The anatomical image to be used as a background. If None, the

728

MNI152 T1 1mm template is used.

729

anat_sform : 4x4 ndarray, optional

730

The affine matrix going from the anatomical image voxel space to

731

MNI space. This parameter is not used when the default

732

anatomical is used, but it is compulsory when using an

733

explicite anatomical image.

734

title : string, optional

735

The title dispayed on the figure.

736

figure_num : integer, optional

737

The number of the matplotlib and Mayavi figures used. If None is

738

given, a new figure is created.

739

mask : 3D ndarray, boolean, optional

740

The brain mask. If None, the mask is computed from the map.

741

auto_sign : boolean, optional

742

If auto_sign is True, the sign of the activation is

743

automaticaly computed: negative activation can thus be

744

plotted.

745

746

Notes

747

-----

748

All the 3D arrays are in numpy convention: (x, y, z)

749

750

Cut coordinates are in Talairach coordinates. Warning: Talairach

751

coordinates are (y, x, z), if (x, y, z) are in voxel-ordering

752

convention.

753

"""

573

754

if do3d:

574

755

if do3d == 'offscreen':

575

756

try:

580

761

plotter = plot_map

581

762

else:

582

763

plotter = plot_map_2d

764

if mask is None:

765

mask = compute_mask(map)

583

766

if vmin is None:

584

767

vmin = np.inf

585

768

pvalue = 0.04

595

778

map *= -1

596

779

if cut_coords is None:

597

780

x, y, z = find_cut_coords(map, activation_threshold=vmin)

598

x, y, z = coord_transform(x, y, z, sform)

599

# XXX: in Talairach the order is (y, x, z)

600

cut_coords = (y, x, z)

781

# XXX: Careful with Voxel/MNI ordering

782

y, x, z = coord_transform(x, y, z, sform)

783

cut_coords = (x, y, z)

601

784

plotter(map, sform, vmin=vmin, cut_coords=cut_coords,

602

785

anat=anat, anat_sform=anat_sform, title=title,

603

786

figure_num=figure_num, mask=mask)

604

787

return vmin, cut_coords

605

788

606

789

607

def plot_niftifile(filename, do3d=False, outputname=None, vmin=None,

790

def plot_niftifile(filename, outputname=None, do3d=False, vmin=None,

608

791

cut_coords=None, anat_filename=None, figure_num=None,

609

792

mask_filename=None, auto_sign=True):

610

""" Given a nifti filename, plot an view of it to a file (png by

793

""" Given a nifti filename, plot a view of it to a file (png by

611

794

default).

795

796

Parameters

797

----------

798

filename : string

799

The name of the Nifti file of the map to be plotted

800

outputname : string, optional

801

The file name of the output file created. By default

802

the name of the input file with a png extension is used.

803

do3d : boolean, optional

804

If do3d is True, a 3D plot is created if Mayavi is installed.

805

vmin : float, optional

806

The lower threshold of the positive activation. This

807

parameter is used to threshold the activation map.

808

cut_coords: 3-tuple of floats, optional

809

The MNI coordinates of the cut to perform, in MNI coordinates

810

and order. If None is given, the cut_coords are automaticaly

811

estimated.

812

anat : string, optional

813

Name of the Nifti image file to be used as a background. If None,

814

the MNI152 T1 1mm template is used.

815

title : string, optional

816

The title dispayed on the figure.

817

figure_num : integer, optional

818

The number of the matplotlib and Mayavi figures used. If None is

819

given, a new figure is created.

820

mask_filename : string, optional

821

Name of the Nifti file to be used as brain mask. If None, the

822

mask is computed from the map.

823

auto_sign : boolean, optional

824

If auto_sign is True, the sign of the activation is

825

automaticaly computed: negative activation can thus be

826

plotted.

827

828

Notes

829

-----

830

831

Cut coordinates are in Talairach coordinates. Warning: Talairach

832

coordinates are (y, x, z), if (x, y, z) are in voxel-ordering

833

convention.

612

834

"""

835

613

836

if outputname is None:

614

837

outputname = os.path.splitext(filename)[0] + '.png'

615

838

if not os.path.exists(filename):

616

839

raise OSError, 'File %s does not exist' % filename

617

840

618

import nifti

619

nim = nifti.NiftiImage(filename)

841

nim = NiftiImage(filename)

620

842

sform = nim.sform

621

843

if any(np.linalg.eigvals(sform)==0):

622

844

raise SformError, "sform affine is not inversible"

623

845

if anat_filename is not None:

624

anat_im = nifti.NiftiImage(anat_filename)

846

anat_im = NiftiImage(anat_filename)

625

847

anat = anat_im.data.T

626

848

anat_sform = anat_im.sform

627

849

else:

629

851

anat_sform = None

630

852

631

853

if mask_filename is not None:

632

mask_im = nifti.NiftiImage(mask_filename)

854

mask_im = NiftiImage(mask_filename)

633

855

mask = mask_im.data.T.astype(np.bool)

634

856

if not np.allclose(mask_im.sform, sform):

635

857

raise SformError, 'Mask does not have same sform as image'

657

879

else:

658

880

fmt = '%s_%04i%s'

659

881

if mask is None:

660

from fff2.utils.mask import compute_mask_intra_array

661

mask = compute_mask_intra_array(nim.data.mean(axis=0)).T

882

mask = compute_mask(nim.data.mean(axis=0)).T

662

883

for index, data in enumerate(nim.data):

663

884

map = data.T

664

885

auto_plot_map(map, sform, vmin=vmin, cut_coords=cut_coords,

678

899

test_coord_transform_trivial()

679

900

test_find_cut_coords()

680

901

681

import sys

682

902

filename = sys.argv[1]

683

903

print "Rendering a visualization of %s" % filename

684

904

plot_niftifile(filename, do3d=True)