Merge lp:~nipy-developers/nipy/trunk-mb into lp:~nipy-developers/nipy/obsolete-do-not-use

[Matthew Brett]

[Matthew Brett (matthew-brett)]

Some updates to io.imageformats from upstream; restored hard dependency checks; fixed missing ellipsis coordmap slicing.

[lp:~nipy-developers/nipy/trunk-mb]

1811. By Matthew Brett <email address hidden> on 2010-01-13

Double ellipsis test

1812. By Matthew Brett <email address hidden> on 2010-01-13

Reworking of Ellipsis code for clarity

[Matthew Brett (matthew-brett)]

I think Gael's approved this one on the list...

review: Approve

[Fernando Perez (fdo.perez)]

Very minor nits

27 +from nipy.core.reference.coordinate_map import CoordinateMap, Affine, compose, \
28 + drop_io_dim, append_io_dim

I'd use

from ... import (...
...)

now that this syntax is allowed.

In

129 +@parametric
130 +def test_array_coord_map():

Maybe a few comments? If one of those tests ever breaks, someone who understands the code less could be in a tight spot to fix them back.

In:

345 +def load(filename, *args, **kwargs):
346 + ''' Load file given filename, guessing at file type
347
348 Parameters
349 ----------
350 - filespec : string or file-like
351 + filename : string or file-like

I actualy found 'filespec' to be clearer, because when I see 'filename' I immediately think "it's already a string", where as filespec to me conveys the idea of name-or-open-handle better.

Looks great otherwise, and these are very minor points. Thanks!

[Matthew Brett (matthew-brett)]

Hi,

Thanks very much for the review, it's very good to have some feedback
on the code.

On Wed, Jan 13, 2010 at 4:51 PM, Fernando Perez <email address hidden> wrote:
> Very minor nits
>
> 27      +from nipy.core.reference.coordinate_map import CoordinateMap, Affine, compose, \
> 28      + drop_io_dim, append_io_dim
>
> I'd use
>
> from ... import (...
> ...)
>
> now that this syntax is allowed.

Yes, we should move that way now, I agree.

> 129     +@parametric
> 130     +def test_array_coord_map():
>
> Maybe a few comments?  If one of those tests ever breaks, someone who understands the code less could be in a tight spot to fix them back.

Yes, fair point.

> In:
>
> 345     +def load(filename, *args, **kwargs):
> 346     + ''' Load file given filename, guessing at file type
> 347
> 348     Parameters
> 349     ----------
> 350     - filespec : string or file-like
> 351     + filename : string or file-like
>
> I actualy found 'filespec' to be clearer, because when I see 'filename' I immediately think "it's already a string", where as filespec to me conveys the idea of name-or-open-handle better.

Ah - yes - actually it's much worse than that. In your Ipython
coding frenzy you may have missed that:

img = Spm99AnalyzeImage(arr, aff)
img.to_filename('test')

actually saves to 'test.img', 'test.hdr', and 'test.mat'. I think
that filespec is better too, but Gael was finding the name a little
confusing, because you usually use it just with a filename.

See you,

Matthew

[lp:~nipy-developers/nipy/trunk-mb]

1813. By Matthew Brett <email address hidden> on 2010-01-13

Changes to basis projections in standardized case, as discussed with Jonathan

1814. By Matthew Brett <email address hidden> on 2010-01-13

Cosmetic code changes to tests_pca

1815. By Matthew Brett <email address hidden> on 2010-01-14

Small cosmetics for PCA routine

[Matthew Brett (matthew-brett)]

Hi,

>> In:
>>
>> 345     +def load(filename, *args, **kwargs):
>> 346     + ''' Load file given filename, guessing at file type
>> 347
>> 348     Parameters
>> 349     ----------
>> 350     - filespec : string or file-like
>> 351     + filename : string or file-like
>>
>> I actualy found 'filespec' to be clearer, because when I see 'filename' I immediately think "it's already a string", where as filespec to me conveys the idea of name-or-open-handle better.
>
> Ah - yes - actually it's much worse than that.   In your Ipython
> coding frenzy you may have missed that:
>
> img = Spm99AnalyzeImage(arr, aff)
> img.to_filename('test')
>
> actually saves to 'test.img', 'test.hdr', and 'test.mat'.   I think
> that filespec is better too, but Gael was finding the name a little
> confusing, because you usually use it just with a filename.

Thinking about this further, I am less sure that the change is bad,
that is, 'to_filename' may indeed be preferable. 'to_filename' does
in fact require a string, because for most formats, more than one file
is needed to save the image, so file object would not work for these
formats. The reason I initially chose 'to_filespec' was to try and
signal the oddness above, but, given that 'to_filespec' made you think
you could use a file object, perhaps it's best left as 'to_filename',
and we'll hope the oddness doesn't show up too much...

See you,

Matthew

[lp:~nipy-developers/nipy/trunk-mb]

1816. By Matthew Brett <email address hidden> on 2010-01-15

PCA standardize with np.std instead of error sum of squares

1817. By Matthew Brett <email address hidden> on 2010-01-15

Maybe more obvious standardize pattern

1818. By Matthew Brett <email address hidden> on 2010-01-16

Cosmetic PCA doc edits

1819. By Matthew Brett <email address hidden> on 2010-01-16

Applied upstream patch continuing to_filespec -> to_filename rename

1820. By Matthew Brett <email address hidden> on 2010-01-16

Moved from deprecated to_filespec method

1821. By Matthew Brett <email address hidden> on 2010-01-16

Removed some remaining references to neuroimaging pacage

1822. By Matthew Brett <email address hidden> on 2010-01-16

Commenting on array_coords test

1823. By Matthew Brett <email address hidden> on 2010-01-16

Removed outdated pynifti docs

1824. By Matthew Brett <email address hidden> on 2010-01-17

Removed old references to pynifti; added affine orientation calculation

[Gael Varoquaux (gael-varoquaux)]

A few remarks:

In pca.py, you are using numpy linalg. I tend to frown on that, because my experience is that quite often numpy linalg is compiled using the included mini-blas, which is not terribly good (for instance, at neurospin, it is the case, but I've seen it elsewhere). As we have a dependency on scipy, if you use scipy.linalg, you know that you are relying on good quality linear algebra libraries.

Also, I see that you use np.std. I remember moving away from it, because it had very bad performance on big arrays:
{{{
In [6]: t = np.random.random((50000, 800))

In [7]: %timeit np.sqrt((np.square(t).sum(axis=0))/t.shape[0])
1 loops, best of 3: 1.19 s per loop

In [8]: %timeit np.sqrt((np.square(t).sum(axis=1))/t.shape[0])
1 loops, best of 3: 587 ms per loop

In [9]: %timeit t.std(axis=1)
1 loops, best of 3: 1.19 s per loop

In [10]: %timeit t.std(axis=0)
1 loops, best of 3: 2.33 s per loop
}}}

[Gael Varoquaux (gael-varoquaux)]

I just remembered, about std, the killer was not the speed, but the memory usage:

{{{
In [1]: import numpy as np

In [2]: t = np.random.random((50000, 3000))

In [3]: %timeit t.std(axis=0)

---------------------------------------------------------------------------
MemoryError Traceback (most recent call last)
}}}

{{{
In [1]: import numpy as np

In [2]: t = np.random.random((50000, 3000))

In [3]: %timeit np.sqrt((np.square(t).sum(axis=0))/t.shape[0])
1 loops, best of 3: 5.03 s per loop
}}}

With my data, std simply wouldn't work.

[Matthew Brett (matthew-brett)]

HI,

Thanks - both the linalg and the std experience is useful, I'll change
the code...

Matthew

On Sun, Jan 17, 2010 at 9:31 AM, Gael Varoquaux
<email address hidden> wrote:
> I just remembered, about std, the killer was not the speed, but the memory usage:
>
> {{{
> In [1]: import numpy as np
>
> In [2]: t = np.random.random((50000, 3000))
>
> In [3]: %timeit t.std(axis=0)
>
> ---------------------------------------------------------------------------
> MemoryError                               Traceback (most recent call last)
> }}}
>
> {{{
> In [1]: import numpy as np
>
> In [2]: t = np.random.random((50000, 3000))
>
> In [3]: %timeit np.sqrt((np.square(t).sum(axis=0))/t.shape[0])
> 1 loops, best of 3: 5.03 s per loop
> }}}
>
> With my data, std simply wouldn't work.
> --
> https://code.launchpad.net/~nipy-developers/nipy/trunk-mb/+merge/17287
> You proposed lp:~nipy-developers/nipy/trunk-mb for merging.
>

[Matthew Brett (matthew-brett)]

Ah, except of course std(t, axis=0) is not the same as
np.sqrt((np.square(t).sum(axis=0))/t.shape[0]) unless np.all(mean(t)
== 0) - which might explain the speed difference.

On Sun, Jan 17, 2010 at 9:31 AM, Gael Varoquaux
<email address hidden> wrote:
> I just remembered, about std, the killer was not the speed, but the memory usage:
>
> {{{
> In [1]: import numpy as np
>
> In [2]: t = np.random.random((50000, 3000))
>
> In [3]: %timeit t.std(axis=0)
>
> ---------------------------------------------------------------------------
> MemoryError                               Traceback (most recent call last)
> }}}
>
> {{{
> In [1]: import numpy as np
>
> In [2]: t = np.random.random((50000, 3000))
>
> In [3]: %timeit np.sqrt((np.square(t).sum(axis=0))/t.shape[0])
> 1 loops, best of 3: 5.03 s per loop
> }}}
>
> With my data, std simply wouldn't work.
> --
> https://code.launchpad.net/~nipy-developers/nipy/trunk-mb/+merge/17287
> You proposed lp:~nipy-developers/nipy/trunk-mb for merging.
>

[lp:~nipy-developers/nipy/trunk-mb]

1825. By Matthew Brett <email address hidden> on 2010-01-17

Use scipy.linalg; do use root MSE instead of standard deviation to standardize - thanks to GV

[Fernando Perez (fdo.perez)]

Looks fine from here, just one question: if we now use 'filename' everywhere, do we still accept open file handles? If so, then we should say somewhere clearly (probably in every docstring) something like:

filename : string or open file-like object

If we accept name-or-object in a variable called 'name', this should be made clear at least in the docstrings (since we've chosen to use deliberately misleading naming convention, we should at least give users a chance of understand the trick and undo the confusion in their heads via documentation)

[Matthew Brett (matthew-brett)]

Hi,

On Mon, Jan 18, 2010 at 1:45 AM, Fernando Perez <email address hidden> wrote:
> Looks fine from here, just one question: if we now use 'filename' everywhere, do we still accept open file handles?  If so, then we should say somewhere clearly (probably in every docstring) something like:
>
>
> filename : string or open file-like object

At the moment, there's a different routine for writing to filenames
and / or file objects - called - perhaps confusingly - img.to_files

The reason there's a different routine is that several image formats
need more than one file-object to write to. For example, you have to
write SPM analyze images by writing filename.img for the binary data
and filename.hdr for the header. If you wanted to do that, passing a
single file object can't work. I could see other ways round this -
like - raising an error with a single file-handle for img.to_filename
when the format needs two files, but it would make the interface a bit
difficult to explain - at least - that was my thought...

See you,

Matthew

[Fernando Perez (fdo.perez)]

On Sun, Jan 17, 2010 at 10:33 PM, Matthew Brett <email address hidden> wrote:
> The reason there's a different routine is that several image formats
> need more than one file-object to write to.  For example, you have to
> write SPM analyze images by writing filename.img for the binary data
> and filename.hdr for the header.   If you wanted to do that, passing a
> single file object can't work.     I could see other ways round this -
> like - raising an error with a single file-handle for img.to_filename
> when the format needs two files, but it would make the interface a bit
> difficult to explain - at least - that was my thought...

OK, then it makes sense to allow *only* actual filenames, and never
open handles on these functions, I see.

I don't think this is worth spending much more time on, as long as
this behavior is explained in the docstrings, since there's a
perfectly valid reason for it, we're good.

Cheers,

f

[lp:~nipy-developers/nipy/trunk-mb]

1826. By Matthew Brett <email address hidden> on 2010-01-18

Applied upstream patch completing image reorientation

1827. By Matthew Brett <email address hidden> on 2010-01-19

Jonathans generalization of affine orientation

1828. By Matthew Brett <email address hidden> on 2010-01-19

Upstream rejigging of orientations

1829. By Matthew Brett <email address hidden> on 2010-01-21

Fix reordering - I was doing it the wrong way round

1830. By Matthew Brett <email address hidden> on 2010-01-21

Montages draft

1831. By Matthew Brett <email address hidden> on 2010-01-21

Rot90 for more obvious slice display

1832. By Matthew Brett <email address hidden> on 2010-01-21

Yoked viz

1833. By Matthew Brett <email address hidden> on 2010-01-21

Reconsidered viz; a bit too drafty, needs proper thought

[Matthew Brett (matthew-brett)]

Any more reviews here? I'll merge tomorrow unless I hear otherwise...

review: Needs Information

=== modified file 'doc/devel/guidelines/coverage_testing.rst'
--- doc/devel/guidelines/coverage_testing.rst	2009-03-09 19:13:28 +0000
+++ doc/devel/guidelines/coverage_testing.rst	2010-01-21 22:28:12 +0000
@@ -28,15 +28,15 @@
 
  Name                                            Stmts   Exec  Cover   Missing
  -----------------------------------------------------------------------------
- neuroimaging                                       21     14    66%   70-74, 88-89
- neuroimaging.core                                   4      4   100%   
- neuroimaging.core.reference                         8      8   100%   
- neuroimaging.core.reference.array_coords          100     90    90%   133-134, 148-151, 220, 222, 235, 242
- neuroimaging.core.reference.coordinate_map        188    187    99%   738
- neuroimaging.core.reference.coordinate_system      61     61   100%   
- neuroimaging.core.reference.slices                 34     34   100%   
- neuroimaging.core.transforms                        0      0   100%   
- neuroimaging.core.transforms.affines               14     14   100%   
+ nipy                                       21     14    66%   70-74, 88-89
+ nipy.core                                   4      4   100%   
+ nipy.core.reference                         8      8   100%   
+ nipy.core.reference.array_coords          100     90    90%   133-134, 148-151, 220, 222, 235, 242
+ nipy.core.reference.coordinate_map        188    187    99%   738
+ nipy.core.reference.coordinate_system      61     61   100%   
+ nipy.core.reference.slices                 34     34   100%   
+ nipy.core.transforms                        0      0   100%   
+ nipy.core.transforms.affines               14     14   100%   
 
 
 The coverage report will cover any python source module imported after
@@ -48,7 +48,7 @@
 ``--cover-package``.  For example, in writing tests for the
 coordinate_map module::
 
-    nosetests --with-coverage --cover-package=neuroimaging.core.reference.coordinate_map test_coordinate_map.py
+    nosetests --with-coverage --cover-package=nipy.core.reference.coordinate_map test_coordinate_map.py
 
 Since that's a lot to type, I wrote a tool called ``sneeze`` to that
 simplifies coverage testing with nose.

=== modified file 'doc/devel/index.rst'
--- doc/devel/index.rst	2009-02-13 02:17:31 +0000
+++ doc/devel/index.rst	2010-01-21 22:28:12 +0000
@@ -19,5 +19,3 @@
    software_design/index
    software_design/usecases/index
    tools/index
-
-   pynifti

=== modified file 'doc/devel/install/ubuntu.rst'
--- doc/devel/install/ubuntu.rst	2009-08-05 02:03:11 +0000
+++ doc/devel/install/ubuntu.rst	2010-01-21 22:28:12 +0000
@@ -13,12 +13,12 @@
    sudo apt-get install python-numpy python-numpy-dev python-scipy
    sudo apt-get install liblapack-dev
    sudo apt-get install python-sympy
-   sudo apt-get install python-nifti
 
 Options::
 
    sudo apt-get install ipython
    sudo apt-get install python-matplotlib
+   sudo apt-get install mayavi2
 
 For getting the code via version control::
 

=== removed file 'doc/devel/pynifti.rst'
--- doc/devel/pynifti.rst	2009-02-03 08:08:16 +0000
+++ doc/devel/pynifti.rst	1970-01-01 00:00:00 +0000
@@ -1,167 +0,0 @@
-.. _pynifti-io:
-
-============
- Pynifti IO
-============
-
-We are using pynifti_ for the underlying nifti file I/O support.
-Pynifti is built upon the nifticlibs_ so this will provide us will
-*full* nifti support.  We are working closely with the author of
-pynifti_, Michael Hanke, and pushing any bug fixes or feature
-improvements upstream to the git repository.
-
-Developers should read through Michael's documentation on the pynifti_
-site for some details on the project.  The source is checked out from
-the `git repository.
-<http://git.debian.org/?p=pkg-exppsy/pynifti.git>`_
-
-Using the command::
-
-  git clone http://git.debian.org/git/pkg-exppsy/pynifti.git
-
-Git
----
-
-Pynifti uses git_ for it's version control system.  Git is very
-different from `svn <http://subversion.tigris.org/>`_, developers
-should read some documentation on git before doing any work with the
-git repository.
-
-The git_ website has several tutorials and full documentation.  A good
-starting point may be the `Git for SVN Users
-<http://git.or.cz/course/svn.html>`_
-
-Git has a unique mechanism for storing multiple branches on your
-machine.  Instead of having separate file directories, git will store
-all branches in one directory and store *branch diffs* in an internal
-database.  When you switch branches (``git checkout branch-name``),
-git will apply the branch diff to the directory, updating any files to
-match the new branch.
-
-Git uses man pages for it's installed documentation.  As with all man
-pages, these contain a lot of useful information, so you should know
-how to access them.  All git commands can be called in two forms:
-
-1. git add <filename>
-
-2. git-add <filename>
-
-The first form is the one you will probably use most and is what is
-often shown in the documentation.  The second form, however, is what
-you need to access the man page.  To see the man page on how to add a
-file to a git repository::
-
-  man git-add
-
-To see a list of all git commands look at the main git man page::
-
-  man git
-
-As with Bazaar, you should identify yourself to git so the repository
-keeps track of who made your edits::
-
-  git config --global user.name "Your Name Comes Here"
-  git config --global user.email [e-mail]
-
-To list your git configuration::
-
-    cburns@pynifti 13:32:31 $ git config -l
-    user.name=Christopher Burns
-    user.email=cburns[at]berkeley[dot]edu
-    color.diff=auto
-    color.status=auto
-    core.repositoryformatversion=0
-    core.filemode=true
-    core.bare=false
-    core.logallrefupdates=true
-    remote.origin.url=ssh://git.debian.org/git/pkg-exppsy/pynifti.git
-    remote.origin.fetch=+refs/heads/*:refs/remotes/origin/*
-    branch.master.remote=origin
-    branch.master.merge=refs/heads/master
-
-We can also see the remote origin branch from which we have cloned.
-
-Development Cycle
------------------
-
-There are 3 development branches that the nipy developers need to
-interact with in the git repository:
-
-* master - the main pynifti repository
-
-* cb/master - Chris' pynifti developer repository
-
-* cb/nipy - Chris' pynifti developer repository with nipy specific code
-
-A nipy development path would look like this:
-
-#. In the master branch, merge with the server to get any updates from
-   other pynifti developers.
-
-#. Checkout the cb/master branch and merge from the local master branch.
-
-#. Make any code edits that should be pushed upstream in the cb/master
-   branch.  Michael cherry-picks changes into the master branch.
-
-#. Checkout the cb/nipy branch and merge from the cb/master.  The
-   cb/nipy branch is used as the source for nipy.
-
-Example:
-^^^^^^^^
-
-I'm in my pynifti source directory::
-
-  cburns@pynifti 13:20:35 $ pwd
-  /Users/cburns/src/pynifti
-
-Use the ``git branch`` command without arguments to see all of your
-local branches.  Below we can see that I'm in my ``cb/master``
-branch::
-  
-  cburns@pynifti 13:20:39 $ git branch
-  * cb/master
-    cb/nipy
-    master
-
-I want to switch to the ``master`` branch and update it with the
-server::
-
-  cburns@pynifti 13:26:08 $ git checkout master
-  Switched to branch "master"
-  cburns@pynifti 13:26:16 $ git branch
-    cb/master
-    cb/nipy
-  * master
-
-Pull from the server to update our master branch::
-
-  cburns@pynifti 13:35:52 $ git pull
-  Password: 
-  Already up-to-date.
-
-Switch into ``cb/master`` and merge with the ``master`` branch.
-Remember, these are not separate directories, git *knows* about the
-other branch by name, so we do not specify a path, we specify a branch
-name.::
-
-  cburns@pynifti 13:36:18 $ git checkout cb/master
-  Switched to branch "cb/master"
-
-  cburns@pynifti 13:38:38 $ git merge master
-  Already up-to-date.
-
-To update the nipy source:
-^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-#. Change to the pynifti directory in the nipy developer trunk::
-
-    cd trunk-dev/neuroimaging/externals/pynifti/utils
-
-#. Run the ``update_source.py`` script to update the source.  This
-   assumes a directory structure the pynifti sources are in the
-   directory ``$HOME/src/pynifti``.  Run the script::
-
-    python update_source.py
-
-
-.. include:: ../links_names.txt

=== modified file 'doc/faq/documentation_faq.rst'
--- doc/faq/documentation_faq.rst	2009-03-19 23:17:13 +0000
+++ doc/faq/documentation_faq.rst	2010-01-21 22:28:12 +0000
@@ -25,8 +25,8 @@
 install.  The error may look something like::
 
   **writing output...** about api/generated/gen
-    api/generated/neuroimaging
-    api/generated/neuroimaging.algorithms.fwhm Format: "png" not
+    api/generated/nipy
+    api/generated/nipy.algorithms.fwhm Format: "png" not
     recognized. Use one of: canon cmap cmapx cmapx_np dia dot eps fig
     hpgl imap imap_np ismap mif mp pcl pic plain plain-ext ps ps2 svg
     svgz tk vml vmlz vtx xdot

=== modified file 'doc/users/coordinate_map.rst'
--- doc/users/coordinate_map.rst	2009-07-03 07:17:10 +0000
+++ doc/users/coordinate_map.rst	2010-01-21 22:28:12 +0000
@@ -22,7 +22,7 @@
 
 
 For more on Coordinate Systems and thier properties
-:mod:`neuroimaging.core.reference.coordinate_system`
+:mod:`nipy.core.reference.coordinate_system`
 
 You can introspect a coordinate map
 

=== modified file 'doc/users/installation.rst'
--- doc/users/installation.rst	2009-10-15 18:28:57 +0000
+++ doc/users/installation.rst	2010-01-21 22:28:12 +0000
@@ -18,16 +18,13 @@
 Must Have
 ^^^^^^^^^
 
-  Python_ 2.4 or later
+  Python_ 2.5 or later
   
   NumPy_ 1.2 or later
 
   SciPy_ 0.7 or later
     Numpy and Scipy are high-level, optimized scientific computing libraries.
 
-  PyNifti_
-    We are using pynifti for the underlying file IO for nifti files.
-
   gcc_
     NIPY does contain a few C extensions for optimized
     routines. Therefore, you must have a compiler to build from

=== modified file 'doc/www/conf.py'
--- doc/www/conf.py	2009-02-12 23:26:12 +0000
+++ doc/www/conf.py	2010-01-21 22:28:12 +0000
@@ -28,15 +28,15 @@
 # The objects.inv file has this info for each mapped object:
 #   label-name classifier path-to-html
 # Examples:
-#   neuroimaging.core.image.image.Image class api/generated/neuroimaging.core.image.image.html
-#   neuroimaging.core.image.generators mod api/generated/neuroimaging.core.image.generators.html
+#   nipy.core.image.image.Image class api/generated/nipy.core.image.image.html
+#   nipy.core.image.generators mod api/generated/nipy.core.image.generators.html
 #
 #intersphinx_mapping = {'../html/doc/manual/html': '../build/html/objects.inv'}
 #
 # In reST documents, I can then link to Python objects in the API like this:
 #
-#This is the image class: :class:`neuroimaging.core.image.image`
-#This is the Image.affine method: :meth:`neuroimaging.core.image.image.Image.affine`
+#This is the image class: :class:`nipy.core.image.image`
+#This is the Image.affine method: :meth:`nipy.core.image.image.Image.affine`
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

=== modified file 'nipy/core/api.py'
--- nipy/core/api.py	2009-12-17 20:42:37 +0000
+++ nipy/core/api.py	2010-01-21 22:28:12 +0000
@@ -1,16 +1,24 @@
 """
-Pseudo-package for all of the core symbols from the image object and its reference
-system.  Use this module for importing core names into your namespace. For example:
- from neuorimaging.core.api import Image
+Pseudo-package for all of the core symbols from the image object and its
+reference system.  Use this module for importing core names into your
+namespace.
+
+For example:
+
+>>> from nipy.core.api import Image
 """
 
 # Note: The order of imports is important here.
-from nipy.core.reference.coordinate_system import CoordinateSystem
-from nipy.core.reference.coordinate_map import CoordinateMap, Affine, compose
-from nipy.core.reference.array_coords import Grid, ArrayCoordMap
-
-from nipy.core.image.image import Image, merge_images, fromarray, is_image
-
-from nipy.core.image.image_list import ImageList
-
-from nipy.core.image.generators import parcels, data_generator, write_data, slice_generator, f_generator, matrix_generator
+from .reference.coordinate_system import CoordinateSystem
+from .reference.coordinate_map import (CoordinateMap, Affine, compose, 
+                                       drop_io_dim, append_io_dim)
+from .reference.array_coords import Grid, ArrayCoordMap
+
+from .image.image import Image, merge_images, fromarray, is_image
+
+from .image.image_list import ImageList
+
+from .image.generators import (parcels, data_generator, write_data,
+                               slice_generator, f_generator,
+                               matrix_generator)
+

=== modified file 'nipy/core/image/image.py'
--- nipy/core/image/image.py	2009-12-17 20:42:37 +0000
+++ nipy/core/image/image.py	2010-01-21 22:28:12 +0000
@@ -76,7 +76,7 @@
         # ii) the shapes are consistent
 
         # self._data is an array-like object.  It must implement a subset of
-        # array methods  (Need to specify these, for now implied in pyniftio)
+        # array methods  (Need to specify these)
         self._data = data
         self._coordmap = coordmap
 

=== modified file 'nipy/core/reference/array_coords.py'
--- nipy/core/reference/array_coords.py	2009-02-25 04:54:32 +0000
+++ nipy/core/reference/array_coords.py	2010-01-21 22:28:12 +0000
@@ -57,11 +57,9 @@
         values : array
            Values of self.coordmap evaluated at np.indices(self.shape).
         """
-
         indices = np.indices(self.shape).astype(
             self.coordmap.input_coords.coord_dtype)
         tmp_shape = indices.shape
-
         # reshape indices to be a sequence of coordinates
         indices.shape = (self.coordmap.ndim[0], np.product(self.shape))
         _range = self.coordmap(indices.T)
@@ -87,13 +85,27 @@
 
         Parameters
         ----------
-        index : ``int`` or ``slice``
-            sequence of integers or slices
-        
+        index : int or tuple
+           int, or sequence of any combination of integers, slices.  The
+           sequence can also cantian one Ellipsis. 
         """
-
+        # index can be single int or tuple
         if type(index) != type(()):
             index = (index,)
+        # allow slicing of form [...,1]
+        if Ellipsis in index:  
+            if sum([i == Ellipsis for i in index]) > 1:
+                raise ValueError("only one Ellipsis (...) allowed in slice")
+            # convert ellipsis to series of slice(None) objects.  For
+            # example, if the coordmap is length 3, we convert (...,1)
+            # to (slice(None), slice(None), 1) - equivalent to [:,:,1]
+            ellipsis_start = list(index).index(Ellipsis)
+            inds_after_ellipsis = index[(ellipsis_start+1):]
+            # the ellipsis continues until any remaining slice specification
+            n_ellipses = len(self.shape) - ellipsis_start - len(inds_after_ellipsis)
+            index = (index[:ellipsis_start]
+                     + n_ellipses * (slice(None),)
+                     + inds_after_ellipsis)
         return _slice(self.coordmap, self.shape, *index)
 
     @staticmethod
@@ -106,16 +118,15 @@
         slices = tuple([slice(0,s,1) for s in shape])
         return Grid(coordmap)[slices]
 
+
 def _slice(coordmap, shape, *slices):
     """
     Slice a 'voxel' CoordinateMap's input_coords with slices. A
     'voxel' CoordinateMap is interpreted as a coordmap having a shape.
     """
-    
     if len(slices) < coordmap.ndim[0]:
         slices = (list(slices) +
                   [slice(None,None,None)] * (coordmap.ndim[0] - len(slices)))
-
     ranges = [np.arange(s) for s in shape]
     cmaps = []
     keep_in_output = []
@@ -177,6 +188,7 @@
     slice_cmap = Affine(A, input_coords, coordmap.input_coords)
     return ArrayCoordMap(compose(coordmap, slice_cmap), tuple(newshape))
                    
+
 class Grid(object):
     """
     Simple class to construct Affine instances with slice notation

=== added file 'nipy/core/reference/tests/test_array_coords.py'
--- nipy/core/reference/tests/test_array_coords.py	1970-01-01 00:00:00 +0000
+++ nipy/core/reference/tests/test_array_coords.py	2010-01-21 22:28:12 +0000
@@ -0,0 +1,88 @@
+""" Testing array coords
+
+"""
+
+import numpy as np
+
+from nose.tools import assert_true, assert_false, \
+     assert_equal, assert_raises
+
+from numpy.testing import assert_array_equal, assert_array_almost_equal
+
+from nipy.testing import parametric
+
+from nipy.core.api import Affine
+
+import nipy.core.reference.array_coords as acs
+
+
+@parametric
+def test_array_coord_map():
+    # array coord map recreates the affine when you slice an image.  In
+    # general, if you take an integer slice in some dimension, the
+    # corresponding column of the affine will go, leaving a row for the
+    # lost dimension, with all zeros, execpt for the translation in the
+    # now-removed dimension, encoding the position of that particular
+    # slice
+    xz = 1.1; yz = 2.3; zz = 3.5
+    xt = 10.0; yt = 11; zt = 12
+    aff = np.diag([xz, yz, zz, 1])
+    aff[:3,3] = [xt, yt, zt]
+    shape = (2,3,4)
+    cmap = Affine.from_params('ijk', 'xyz', aff)
+    acm = acs.ArrayCoordMap(cmap, shape)
+    # slice the coordinate map for the first axis
+    sacm = acm[1]
+    # The affine has lost the first column, but has a remaining row (the
+    # first) encoding the translation to get to this slice
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [0, 0, xz+xt],
+                [yz, 0, yt],
+                [0, zz, zt],
+                [0, 0, 1]]))
+    sacm = acm[:,1]
+    # lost second column, remaining second row with translation
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [xz, 0, xt],
+                [0, 0, yz+yt],
+                [0, zz, zt],
+                [0, 0, 1]]))
+    sacm = acm[:,:,2]
+    # ditto third column and row
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [xz, 0, xt],
+                [0, yz, yt],
+                [0, 0, 2*zz+zt],
+                [0, 0, 1]]))
+    # check ellipsis slicing is the same as [:,: ...
+    sacm = acm[...,2]
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [xz, 0, xt],
+                [0, yz, yt],
+                [0, 0, 2*zz+zt],
+                [0, 0, 1]]))
+    # that ellipsis can follow other slice types
+    sacm = acm[:,...,2]
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [xz, 0, xt],
+                [0, yz, yt],
+                [0, 0, 2*zz+zt],
+                [0, 0, 1]]))
+    # that there can be only one ellipsis
+    yield assert_raises(ValueError, acm.__getitem__, (
+            (Ellipsis,Ellipsis,2)))
+    # that you can integer slice in all three dimensions, leaving only
+    # the translation column
+    sacm = acm[1,0,2]
+    yield assert_array_almost_equal(sacm.coordmap.affine,
+                             np.array([
+                [xz+xt],
+                [yt],
+                [2*zz+zt],
+                [1]]))
+    

=== modified file 'nipy/io/files.py'
--- nipy/io/files.py	2009-12-17 20:42:37 +0000
+++ nipy/io/files.py	2010-01-21 22:28:12 +0000
@@ -207,7 +207,7 @@
     # Set zooms
     hdr.set_zooms(zooms)
     # save to disk
-    out_img.to_filespec(filename)
+    out_img.to_filename(filename)
     return Fimg
 
 

=== modified file 'nipy/io/imageformats/__init__.py'
--- nipy/io/imageformats/__init__.py	2009-12-16 22:40:01 +0000
+++ nipy/io/imageformats/__init__.py	2010-01-21 22:28:12 +0000
@@ -2,20 +2,20 @@
 #ex: set sts=4 ts=4 sw=4 et:
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 #
-#   See COPYING file distributed along with the PyNIfTI package for the
+#   See COPYING file distributed along with the NiBabel package for the
 #   copyright and license terms.
 #
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 """This module provides Python bindings to the NIfTI data format.
 
-The PyNIfTI module is a Python interface to the NIfTI I/O libraries. Using
-PyNIfTI, one can easily read and write NIfTI and ANALYZE images from within
+The NiBabel module is a Python interface to the NIfTI I/O libraries. Using
+NiBabel, one can easily read and write NIfTI and ANALYZE images from within
 Python. The :class:`~nifti.image.NiftiImage` class provides pythonic
 access to the full header information and for a maximum of interoperability the
 image data is made available via NumPy arrays.
 
 ===============================
- pynifti python implementation
+ nibabel python implementation
 ===============================
 
 Quickstart::
@@ -59,4 +59,8 @@
 from nipy.io.imageformats.spm2analyze import Spm2AnalyzeHeader, Spm2AnalyzeImage
 from nipy.io.imageformats.nifti1 import Nifti1Header, Nifti1Image
 from nipy.io.imageformats.minc import MincHeader, MincImage
-from nipy.io.imageformats.funcs import squeeze_image, concat_images, four_to_three
+from nipy.io.imageformats.funcs import (squeeze_image, concat_images, four_to_three,
+                                        as_closest_canonical)
+from nipy.io.imageformats.orientations import (io_orientation, orientation_affine,
+                                   flip_axis, OrientationError,
+                                   apply_orientation)

=== modified file 'nipy/io/imageformats/analyze.py'
--- nipy/io/imageformats/analyze.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/analyze.py	2010-01-21 22:28:12 +0000
@@ -19,7 +19,8 @@
 therefore must know how to predict image filenames from header
 filenames, whether these are different, and so on.
 
-You can access and set fields of a particular header type using standard __getitem__ / __setitem__ syntax:
+You can access and set fields of a particular header type using standard
+__getitem__ / __setitem__ syntax:
 
     hdr['field'] = 10
 
@@ -94,13 +95,16 @@
 
 '''
 
+import warnings
+
 import numpy as np
 
 from nipy.io.imageformats.volumeutils import pretty_mapping, endian_codes, \
      native_code, swapped_code, hdr_getterfunc, \
      make_dt_codes, HeaderDataError, HeaderTypeError, allopen
 
-from nipy.io.imageformats.header_ufuncs import read_data, write_data, adapt_header
+from nipy.io.imageformats.header_ufuncs import read_data, write_data, \
+    adapt_header, read_unscaled_data
 
 from nipy.io.imageformats import imageglobals as imageglobals
 from nipy.io.imageformats.spatialimages import SpatialImage
@@ -1071,6 +1075,41 @@
         self._data = read_data(self._header, allopen(fname))
         return self._data
 
+    def get_unscaled_data(self):
+        """ Return image data without image scaling applied
+
+        Summary: please use the ``get_data`` method instead of this
+        method unless you are sure what you are doing, and that you will
+        only be using image formats for which this method exists and
+        returns sensible results.
+        
+        Use this method with care; the modified Analyze-type formats
+        such as SPM formats, and nifti1, specify that the image data
+        array, as they are expecting to return it, is given by the raw
+        data on disk, multiplied by a scalefactor and maybe with the
+        addition of a constant.  This method returns the data on the
+        disk, without these format-specific scalings applied.  Please
+        use this method only if you absolutely need the unscaled data,
+        and the magnitude of the data, as given by the scalefactor, is
+        not relevant to your application.  The Analyze-type formats have
+        a single scalefactor +/- offset per image on disk. If you do not
+        care about the absolute values, and will be removing the mean
+        from the data, then the unscaled values will have preserved
+        intensity ratios compared to the mean-centered scaled data.
+        However, this is not necessarily true of other formats with more
+        complicated scaling - such as MINC.
+
+        Note that - unlike the scaled ``get_data`` method, we do not
+        cache the array, to minimize the memory taken by the object.
+        """
+        if not self._files:
+            return None
+        try:
+            fname = self._files['image']
+        except KeyError:
+            return None
+        return read_unscaled_data(self._header, allopen(fname))
+        
     def get_header(self):
         ''' Return header
 
@@ -1099,11 +1138,6 @@
         self._header.set_data_dtype(dtype)
     
     @classmethod
-    def from_filespec(klass, filespec):
-        files = klass.filespec_to_files(filespec)
-        return klass.from_files(files)
-    
-    @classmethod
     def from_files(klass, files):
         fname = files['header']
         header = klass._header_maker.from_fileobj(allopen(fname))
@@ -1112,14 +1146,6 @@
         ret._files = files
         return ret
     
-    @classmethod
-    def from_image(klass, img):
-        orig_hdr = img.get_header()
-        return klass(img.get_data(),
-                     img.get_affine(),
-                     img.get_header(),
-                     img.extra)
-    
     @staticmethod
     def filespec_to_files(filespec):
         ftups = filetuples.FileTuples(
@@ -1133,12 +1159,6 @@
         files = dict(zip(('header', 'image'), ftups.get_filenames()))
         return files
 
-    def to_filespec(self, filespec):
-        ''' Write image to files given by filespec
-        '''
-        files = self.filespec_to_files(filespec)
-        self.to_files(files)
-
     def to_files(self, files=None):
         ''' Write image to files passed, or self._files
         '''
@@ -1196,16 +1216,7 @@
             RZS = self._affine[:3,:3]
             vox = np.sqrt(np.sum(RZS * RZS, axis=0))
             hdr['pixdim'][1:4] = vox
-        
-    @classmethod
-    def load(klass, filespec):
-        return klass.from_filespec(filespec)
-
-    @classmethod
-    def save(klass, img, filespec):
-        img = klass.from_image(img)
-        img.to_filespec(filespec)
 
 
 load = AnalyzeImage.load
-save = AnalyzeImage.save
+save = AnalyzeImage.instance_to_filename

=== modified file 'nipy/io/imageformats/funcs.py'
--- nipy/io/imageformats/funcs.py	2009-12-16 22:40:01 +0000
+++ nipy/io/imageformats/funcs.py	2010-01-21 22:28:12 +0000
@@ -1,6 +1,10 @@
 ''' Processor functions for images '''
 import numpy as np
 
+from .orientations import (io_orientation, orientation_affine, flip_axis,
+                            apply_orientation, OrientationError)
+
+
 def squeeze_image(img):
     ''' Return image, remove axes length 1 at end of image shape
 
@@ -119,3 +123,53 @@
         img3d = image_maker(arr3d, affine, header)
         imgs.append(img3d)
     return imgs
+
+
+def as_closest_canonical(img, enforce_diag=False):
+    ''' Return `img` with data reordered to be closest to canonical
+
+    Canonical order is the ordering of the output axes.
+
+    Parameters
+    ----------
+    img : ``spatialimage``
+    enforce_diag : {False, True}, optional
+       If True, before transforming image, check if the resulting image
+       affine will be close to diagonal, and if not, raise an error
+    
+    Returns
+    -------
+    canonical_img : ``spatialimage``
+       Version of `img` where the underlying array may have been
+       reordered and / or flipped so that axes 0,1,2 are those axes in
+       the input data that are, respectively, closest to the output axis
+       orientation.  We modify the affine accordingly.  If `img` is
+       already has the correct data ordering, we just return `img`
+       unmodified.
+    '''
+    aff = img.get_affine()
+    ornt = io_orientation(aff)
+    if np.all(ornt == [[0,1],
+                       [1,1],
+                       [2,1]]): # canonical already
+        # however, the affine may not be diagonal
+        if enforce_diag and not _aff_is_diag(aff):
+            raise OrientationError('Transformed affine is not diagonal')
+        return img
+    shape = img.get_shape()
+    t_aff = orientation_affine(ornt, shape)
+    out_aff = np.dot(aff, t_aff)
+    # check if we are going to end up with something diagonal
+    if enforce_diag and not _aff_is_diag(aff):
+        raise OrientationError('Transformed affine is not diagonal')
+    # we need to transform the data
+    arr = img.get_data()
+    t_arr = apply_orientation(arr, ornt)
+    return img.__class__(t_arr, out_aff, img.get_header())
+
+
+def _aff_is_diag(aff):
+    ''' Utility function returning True if affine is nearly diagonal '''
+    rzs_aff = aff[:3,:3]
+    return np.allclose(rzs_aff, np.diag(np.diag(rzs_aff)))
+

=== removed file 'nipy/io/imageformats/images.py'
--- nipy/io/imageformats/images.py	2009-06-05 03:13:27 +0000
+++ nipy/io/imageformats/images.py	1970-01-01 00:00:00 +0000
@@ -1,124 +0,0 @@
-#emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
-#ex: set sts=4 ts=4 sw=4 et:
-''' Base class for images
-
-A draft.
-'''
-
-import nipy.io.imageformats.ioimps as ioimps
-
-class ImageError(Exception):
-    pass
-
-class Image(object):
-    ''' Base class for images
-
-    Attributes:
-
-    * data : array-like; image data
-    * affine : array; mapping from image voxels to output space
-      coordinates
-    * output_space : string; name for affine output space
-    * meta : dict-like; image metadata
-    * io : image io implementation object
-
-    Properties (sorry guys):
-
-    * filename : string : read only, filename or None if none
-    * mode : string; 'r'ead, 'w'rite, or 'rw' 
-    
-    Methods:
-
-    * save(filespec=None)
-    * get_filespec()
-    * set_filespec(filespec) 
-    * same_as_filename()
-
-    Class attributes:
-
-    * default_io_class
-
-    '''
-
-    default_io_class = ioimps.default_io
-    
-    def __init__(self, data,
-                 affine=None,
-                 output_space=None,
-                 meta=None,
-                 filespec=None,
-                 io=None,
-                 mode='rw'):
-        self.data = data
-        self.affine = affine
-        self.output_space = output_space
-        if meta is None:
-            meta = {}
-        self.meta = meta
-        if not filespec is None:
-            if io is None:
-                io = ioimps.guessed_imp(filespec)
-            else:
-                io.set_filespec(filespec)
-        if io is None:
-            io = self.default_io_class()
-        self.io = io
-        # lay down flag to show changes in IO
-        self._io_hash = io.get_hash()
-        self._mode = None
-        self.mode = mode
-        
-    @property
-    def filename(self):
-        filespec = self.get_filespec()
-        return filespec.filename
-
-    def mode():
-        def fget(self):
-            return self._mode
-        def fset(self, mode):
-            if not set('rw').issuperset(set(mode)):
-                raise ImageError('Invalid mode "%s"' % mode)
-            self._mode = mode
-        doc = 'image read / write mode'
-        return locals()
-    mode = property(**mode())
-            
-    def get_filespec(self):
-        return self.io.get_filespec()
-
-    def set_filespec(self, filespec):
-        self.io.set_filespec(filespec)
-        
-    def save(self, filespec=None, io=None):
-        io.save_image(self, filespec, io)
-
-    @classmethod
-    def from_io(klass, io):
-        return io.as_image(klass)
-                     
-    def same_as_filename(self, filename=None):
-        ''' True if image in memory is same as image on disk '''
-        fname = self.filename
-        if fname is None:
-            return False
-        if filename:
-            if filename != fname:
-                return False
-        # first, check if io has changed
-        if self._io_hash != self.io.get_hash():
-            return False
-        # Then get io to check image against itself
-        return self.io.same_as_image(self)
-    
-
-def load(filespec, maker=Image, io=None):
-    if io is None:
-        io = ioimps.guessed_implementation(filespec)
-    else:
-        io.set_filespec(filespec)
-    return maker.from_io(io)
-
-
-def save(img, filespec=None, io=None):
-    img.save(filespec, io)

=== removed file 'nipy/io/imageformats/ioimps.py'
--- nipy/io/imageformats/ioimps.py	2009-06-05 03:13:27 +0000
+++ nipy/io/imageformats/ioimps.py	1970-01-01 00:00:00 +0000
@@ -1,64 +0,0 @@
-''' IO implementatations '''
-
-def guessed_imp(filespec):
-    ''' Return implementation guessed from filespec '''
-    raise NotImplementedError
-
-class IOImplementation(object):
-    def __init__(self, filespec = None):
-        self._filespec = None
-        self.set_filespec(filespec)
-        
-    def get_filespec(self):
-        return self._filespec
-
-    def set_filespec(self, filespec):
-        self._filespec = filespec
-
-    def to_filespec(self, filespec=None):
-        raise NotImplementedError
-
-    def copy(self):
-        raise NotImplementedError
-
-    def get_affine(self):
-        raise NotImplementedError
-
-    def set_affine(self, affine):
-        raise NotImplementedError
-
-    def get_output_space(self):
-        raise NotImplementedError
-
-    def set_output_space(self, output_space):
-        raise NotImplementedError
-
-    def get_data_shape(self):
-        raise NotImplementedError
-
-    def set_data_shape(self, shape):
-        raise NotImplementedError
-
-    def get_data_dtype(self):
-        raise NotImplementedError
-
-    def set_data_dtype(self, dtype):
-        raise NotImplementedError
-
-    def write_slice(data, slicedef, outfile = None):
-        raise NotImplementedError
-
-    def as_image(self, image_maker):
-        raise NotImplementedError
-
-    def save_image(self, image, filespec=None, io=None):
-        raise NotImplementedError
-
-    def same_as_image(self, image):
-        raise NotImplementedError
-
-    def get_hash(self):
-        raise NotImplementedError
-
-default_io = IOImplementation
-

=== modified file 'nipy/io/imageformats/loadsave.py'
--- nipy/io/imageformats/loadsave.py	2009-06-05 03:13:27 +0000
+++ nipy/io/imageformats/loadsave.py	2010-01-21 22:28:12 +0000
@@ -5,12 +5,12 @@
 from nipy.io.imageformats import minc
 
 
-def load(filespec, *args, **kwargs):
-    ''' Load file given filespec, guessing at file type
+def load(filename, *args, **kwargs):
+    ''' Load file given filename, guessing at file type
 
     Parameters
     ----------
-    filespec : string or file-like
+    filename : string or file-like
        specification of filename or file to load
     *args
     **kwargs
@@ -23,31 +23,31 @@
 
     '''
     # Try and guess file type from filename
-    if isinstance(filespec, basestring):
-        fname = filespec
+    if isinstance(filename, basestring):
+        fname = filename
         for ending in ('.gz', '.bz2'):
-            if filespec.endswith(ending):
+            if filename.endswith(ending):
                 fname = fname[:-len(ending)]
                 break
         if fname.endswith('.nii'):
-            return nifti1.load(filespec, *args, **kwargs)
+            return nifti1.load(filename, *args, **kwargs)
         if fname.endswith('.mnc'):
-            return minc.load(filespec, *args, **kwargs)
+            return minc.load(filename, *args, **kwargs)
     # Not a string, or not recognized as nii or mnc
     try:
-        files = nifti1.Nifti1Image.filespec_to_files(filespec)
+        files = nifti1.Nifti1Image.filespec_to_files(filename)
     except ValueError:
         raise RuntimeError('Cannot work out file type of "%s"' %
-                           filespec)
+                           filename)
     hdr = nifti1.Nifti1Header.from_fileobj(
         vu.allopen(files['header']),
         check=False)
     magic = hdr['magic']
     if magic in ('ni1', 'n+1'):
-        return nifti1.load(filespec, *args, **kwargs)
-    return spm2.load(filespec, *args, **kwargs)
-
-
-def save(img, filespec):
+        return nifti1.load(filename, *args, **kwargs)
+    return spm2.load(filename, *args, **kwargs)
+
+
+def save(img, filename):
     ''' Save an image to file without changing format'''
-    img.to_filespec(filespec)
+    img.to_filename(filename)

=== modified file 'nipy/io/imageformats/minc.py'
--- nipy/io/imageformats/minc.py	2009-09-11 01:35:49 +0000
+++ nipy/io/imageformats/minc.py	2010-01-21 22:28:12 +0000
@@ -256,11 +256,6 @@
         return self._header.get_data_dtype()
     
     @classmethod
-    def from_filespec(klass, filespec):
-        files = klass.filespec_to_files(filespec)
-        return klass.from_files(files)
-    
-    @classmethod
     def from_files(klass, files):
         fname = files['image']
         header = klass._header_maker.from_fileobj(allopen(fname))
@@ -269,20 +264,9 @@
         ret._files = files
         return ret
     
-    @classmethod
-    def from_image(klass, img):
-        return klass(img.get_data(),
-                     img.get_affine(),
-                     img.get_header(),
-                     img.extra)
-    
     @staticmethod
     def filespec_to_files(filespec):
         return {'image':filespec}
         
-    @classmethod
-    def load(klass, filespec):
-        return klass.from_filespec(filespec)
-
 
 load = MincImage.load

=== modified file 'nipy/io/imageformats/nifti1.py'
--- nipy/io/imageformats/nifti1.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/nifti1.py	2010-01-21 22:28:12 +0000
@@ -82,7 +82,7 @@
     (512, 'uint16', np.uint16),
     (768, 'uint32', np.uint32),
     (1024,'int64', np.int64),
-    (1280, 'int64', np.uint64),
+    (1280, 'uint64', np.uint64),
     (1536, 'float128', _float128t), # Only numpy defined on 64 bit
     (1792, 'complex128', np.complex128),
     (2048, 'complex256', _complex256t), # 64 bit again
@@ -1458,19 +1458,15 @@
         data = self.get_data()
         # Adapt header to possible two<->one file difference
         is_pair = files['header'] != files['image']
-
         hdr = self.get_header().for_file_pair(is_pair)
-
         # if any extensions, figure out necessary vox_offset for extensions to
         # fit
         if self.extra.has_key('extensions') and len(self.extra['extensions']):
             hdr['vox_offset'] = len(hdr.binaryblock) \
                                 + self.extra['extensions'].get_sizeondisk()
-
         slope, inter, mn, mx = adapt_header(hdr, data)
         hdrf = allopen(files['header'], 'wb')
         hdr.write_to(hdrf)
-
         # write all extensions to file
         # assumes that the file ptr is right after the magic string
         if not self.extra.has_key('extensions'):
@@ -1478,13 +1474,11 @@
             hdrf.write(np.array((0,0,0,0), dtype=np.int8).tostring())
         else:
             self.extra['extensions'].write_to(hdrf)
-
-
         if is_pair:
             imgf = allopen(files['image'], 'wb')
         else: # single file for header and image
             imgf = hdrf
-        # streams like bz2 do not allow seeks, even forward.  We
+            # streams like bz2 do not allow seeks, even forward.  We
             # check where to go, and write zeros up until the data part
             # of the file
             offset = hdr.get_data_offset()
@@ -1520,4 +1514,4 @@
 
 
 load = Nifti1Image.load
-save = Nifti1Image.save
+save = Nifti1Image.instance_to_filename

=== added file 'nipy/io/imageformats/orientations.py'
--- nipy/io/imageformats/orientations.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/orientations.py	2010-01-21 22:28:12 +0000
@@ -0,0 +1,254 @@
+''' Utilities for calculating and applying affine orientations '''
+
+import numpy as np
+import numpy.linalg as npl
+
+
+class OrientationError(Exception):
+    pass
+
+
+def io_orientation(affine, tol=None):
+    ''' Orientation of input axes in terms of output axes for `affine`
+    
+    Valid for an affine transformation from ``m`` dimensions to ``n``
+    dimensions (``affine.shape == (n+1, m+1)``). 
+
+    The calculated orientations can be used to transform associated
+    arrays to best match the output orientations. If ``n`` > ``m``, then
+    some of the output axes should be considered dropped in this
+    orientation.
+
+    Parameters
+    ----------
+    affine : (n+1,m+1) ndarray-like
+       Transformation affine from ``m`` inputs to ``n`` outputs.
+       Usually this will be a shape (4,4) matrix, transforming 3 inputs
+       to 3 outputs, but the code also handles the more general case
+     tol : {None, float}, optional
+        threshold below which SVD values of the affine are considered
+        zero. If `tol` is None, and ``S`` is an array with singular
+        values for `affine`, and ``eps`` is the epsilon value for
+        datatype of ``S``, then `tol` set to ``S.max() * eps``.
+
+    Returns
+    -------
+    orientations : (n,2) ndarray
+       one row per input axis, where the first value in each row is the
+       closest corresponding output axis, and the second value is 1 if
+       the input axis is in the same direction as the corresponding
+       output axis, , and -1 if it is in the opposite direction, to the
+       corresponding output axis.  If a row is [np.nan, np.nan], which
+       can happen when n > m, then this row should be considered
+       dropped.
+    '''
+    affine = np.asarray(affine)
+    n, m = affine.shape[0]-1, affine.shape[1]-1
+    # extract the underlying rotation matrix
+    RZS = affine[:n,:m]
+    zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
+    RS = RZS / zooms
+    # Transform below is polar decomposition, returning the closest
+    # shearless matrix R to RS
+    P, S, Qs = npl.svd(RS)
+    # Threshold the singular values to determine the rank.
+    if tol is None:
+        tol = S.max() * np.finfo(S.dtype).eps
+    keep = (S > tol)
+    R = np.dot(P[:,keep], Qs[keep])
+    # the matrix R is such that np.dot(R,R.T) is projection onto the
+    # columns of P[:,keep] and np.dot(R.T,R) is projection onto the rows
+    # of Qs[keep].  R (== np.dot(R, np.eye(m))) gives rotation of the
+    # unit input vectors to output coordinates.  Therefore, the row
+    # index of abs max R[:,N], is the output axis changing most as input
+    # axis N changes.  In case there are ties, we choose the axes
+    # iteratively, removing used axes from consideration as we go
+    ornt = np.ones((n,2), dtype=np.int8)*np.nan
+    for in_ax in range(m):
+        col = R[:,in_ax]
+        if not np.alltrue(np.equal(col, 0)):
+            out_ax = np.argmax(np.abs(col))
+            ornt[in_ax,0] = out_ax
+            assert col[out_ax] != 0
+            if col[out_ax] < 0:
+                ornt[in_ax,1] = -1
+            else:
+                ornt[in_ax,1] = 1
+            # remove the identified axis from further consideration, by
+            # zeroing out the corresponding row in R
+            R[out_ax,:] = 0
+    return ornt
+
+
+def _ornt_to_affine(orientations):
+    ''' Create affine transformation matrix determined by orientations.
+    
+    This transformation will simply flip, transpose, and possibly drop some
+    coordinates.
+
+    Parameters
+    ----------
+    orientations : (n,2) ndarray
+       one row per input axis, where the first value in each row is the
+       closest corresponding output axis, and the second value is 1 if
+       the input axis is in the same direction as the corresponding
+       output axis, and -1 if it is in the opposite direction, to the
+       corresponding output axis.  If a row has first entry np.nan, then
+       this axis is dropped from the output.
+
+    Returns
+    -------
+    affine : (m+1,n+1) ndarray
+       matrix representing flipping / dropping axes. m is equal to the 
+       number of rows of orientations[:,0] that are not np.nan
+    '''
+    ornt = np.asarray(orientations)
+    n = ornt.shape[0]
+    keep = ~np.isnan(ornt[:,1])
+    # These are the input coordinate axes that do have a matching output
+    # column in the orientation.  That is, if the 2nd row is [np.nan,
+    # np.nan] then the orientation indicates that no output axes of an
+    # affine with this orientation matches the 2nd input coordinate
+    # axis.  This would happen if the 2nd row of the affine that
+    # generated ornt was [0,0,0,*]
+    axes_kept = np.arange(n)[keep] 
+    m = keep.sum(0)
+    # the matrix P represents the affine transform impled by ornt. If
+    # all entries of ornt are not np.nan, then P is square otherwise it
+    # has more columns than rows indicating some coordinates were
+    # dropped
+    P = np.zeros((m+1,n+1))
+    P[-1,-1] = 1
+    for idx in range(m):
+        axs, flip = ornt[axes_kept[idx]]
+        P[idx, axs] = flip
+    return P
+
+
+def apply_orientation(arr, ornt):
+    ''' Apply transformations implied by `ornt` to the first
+    n axes of the array `arr`
+
+    Parameters
+    ----------
+    arr : array-like of data with ndim >= n
+    ornt : (n,2) orientation array
+       orientation transform. ``ornt[N,1]` is flip of axis N of the
+       array implied by `shape`, where 1 means no flip and -1 means
+       flip.  For example, if ``N==0`` and ``ornt[0,1] == -1``, and
+       there's an array ``arr`` of shape `shape`, the flip would
+       correspond to the effect of ``np.flipud(arr)``.  ``ornt[:,0]`` is
+       the transpose that needs to be done to the implied array, as in
+       ``arr.transpose(ornt[:,0])``
+
+    Returns
+    -------
+    t_arr : ndarray
+       data array `arr` transformed according to ornt
+    '''
+    t_arr = np.asarray(arr)
+    ornt = np.asarray(ornt)
+    n = ornt.shape[0]
+    if t_arr.ndim < n:
+        raise OrientationError('Data array has fewer dimensions than '
+                               'orientation')
+    # no coordinates can be dropped for applying the orientations
+    if np.any(np.isnan(ornt[:,0])):
+        raise OrientationError('Cannot drop coordinates when '
+                               'applying orientation to data')
+    shape = t_arr.shape
+    # apply ornt transformations
+    for ax, flip in enumerate(ornt[:,1]):
+        if flip == -1:
+            t_arr = flip_axis(t_arr, axis=ax)
+    full_transpose = np.arange(t_arr.ndim)
+    # ornt indicates the transpose that has occurred - we reverse it
+    full_transpose[:n] = np.argsort(ornt[:,0])
+    t_arr = t_arr.transpose(full_transpose)
+    return t_arr
+
+
+def orientation_affine(ornt, shape):
+    ''' Affine transform resulting from transforms implied in `ornt`
+
+    Imagine you have an array ``arr`` of shape `shape`, and you apply the
+    transforms implied by `ornt` (more below), to get ``tarr``.
+    ``tarr`` may have a different shape ``shape_prime``.  This routine
+    returns the affine that will take a array coordinate for ``tarr``
+    and give you the corresponding array coordinate in ``arr``.  
+
+    Parameters
+    ----------
+    ornt : (n,2) ndarray
+       orientation transform. ``ornt[N,1]` is flip of axis N of the
+       array implied by `shape`, where 1 means no flip and -1 means
+       flip.  For example, if ``N==0`` and ``ornt[0,1] == -1``, and
+       there's an array ``arr`` of shape `shape`, the flip would
+       correspond to the effect of ``np.flipud(arr)``.  ``ornt[:,0]`` is
+       the transpose that needs to be done to the implied array, as in
+       ``arr.transpose(ornt[:,0])``
+
+    shape : length n sequence
+       shape of array you may transform with `ornt`
+
+    Returns
+    -------
+    transformed_affine : (n+1,n+1) ndarray
+       An array ``arr`` (shape `shape`) might be transformed according
+       to `ornt`, resulting in a transformed array ``tarr``.
+       `transformed_affine` is the transform that takes you from array
+       coordinates in ``tarr`` to array coordinates in ``arr``.
+    '''
+    ornt = np.asarray(ornt)
+    n = ornt.shape[0]
+    shape = np.array(shape)[:n]
+    # ornt implies a flip, followed by a transpose.   We need the affine
+    # that inverts these.  Thus we need the affine that first undoes the
+    # effect of the transpose, then undoes the effects of the flip.
+    # ornt indicates the transpose that has occurred to get the current
+    # ordering, relative to canonical, so we just use that
+    undo_reorder = np.eye(n+1)[list(ornt[:,0]) + [n],:]
+    undo_flip = np.diag(list(ornt[:,1]) + [1.0])
+    center_trans = -(shape-1) / 2.0
+    undo_flip[:n,n] = (ornt[:,1] * center_trans) - center_trans
+    return np.dot(undo_flip, undo_reorder)
+
+
+def flip_axis(arr, axis=0):
+    ''' Flip contents of `axis` in array `arr`
+
+    ``flip_axis`` is the same transform as ``np.flipud``, but for any
+    axis.  For example ``flip_axis(arr, axis=0)`` is the same transform
+    as ``np.flipud(arr)``, and ``flip_axis(arr, axis=1)`` is the same
+    transform as ``np.fliplr(arr)``
+
+    Parameters
+    ----------
+    arr : array-like
+    axis : int, optional
+       axis to flip.  Default `axis` == 0
+
+    Returns
+    -------
+    farr : array
+       Array with axis `axis` flipped
+
+    Examples
+    --------
+    >>> a = np.arange(6).reshape((2,3))
+    >>> a
+    array([[0, 1, 2],
+           [3, 4, 5]])
+    >>> flip_axis(a, axis=0)
+    array([[3, 4, 5],
+           [0, 1, 2]])
+    >>> flip_axis(a, axis=1)
+    array([[2, 1, 0],
+           [5, 4, 3]])
+    '''
+    arr = np.asanyarray(arr)
+    arr = arr.swapaxes(0, axis)
+    arr = np.flipud(arr)
+    return arr.swapaxes(axis, 0)
+
+

=== modified file 'nipy/io/imageformats/spatialimages.py'
--- nipy/io/imageformats/spatialimages.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/spatialimages.py	2010-01-21 22:28:12 +0000
@@ -6,57 +6,79 @@
 that is specific to the image format - and ``extra`` - a dictionary
 container for any other metadata.
 
-It has attributes::
+It has attributes:
 
-    extra
-    filename (read only)
+   * extra
     
-and methods::
-
-    .get_data()
-    .get_affine()
-    .get_header()
-    .to_files() # writes image out to passed or
-    .get_raw_data()
-    .write_data(fileobj)
-    .write_raw_data(fileobj)
+methods:
+
+   * .get_data()
+   * .get_affine()
+   * .get_header()
+   * .get_shape()
+   * .set_shape(shape)
+   * .to_filename(fname) - writes data to filename(s) derived from
+     ``fname``, where the derivation may differ between formats.
+   * to_files() - save image to files with which the image is already
+     associated.  Or ``img.to_files(files)`` saves to the files passed.
+
+classmethods:
+
+   * from_filename(fname) - make instance by loading from filename
+   * instance_to_filename(img, fname) - save ``img`` instance to
+     filename ``fname``.
 
 There are several ways of writing data.
 =======================================
 
 There is the usual way, which is the default::
 
-    img.write_data(data, fileobj)
+    img.to_filename(fname)
 
-and that is, to take the data array, ``data``, and cast it to the
-datatype the header expects, setting any available header scaling
+and that is, to take the data encapsulated by the image and cast it to
+the datatype the header expects, setting any available header scaling
 into the header to help the data match.
 
+You can load the data into an image from file with::
+
+   img.from_filename(fname)
+
+The image stores its associated files in a rather secretive way.  In
+order to just save an image, for which you know there is an associated
+filename, or other storage, you can do::
+
+   img.to_files()
+
+alternatively, you can pass in the needed files yourself, into this
+method, as an argument.
+
 You can get the data out again with of::
 
     img.get_data(fileobj)
 
-Less commonly, you might want to fetch out the unscaled array via
-the header::
-
-    unscaled_data = img.get_raw_data(fileobj)
-
-then do something with it.  Then put it back again::
-
-    img.write_raw_data(modifed_unscaled_data, fileobj)
+Less commonly, for some image types that support it, you might want to
+fetch out the unscaled array via the header::
+
+    unscaled_data = img.get_unscaled_data()
+
+Analyze-type images (including nifti) support this, but others may not
+(MINC, for example).
 
 Sometimes you might to avoid any loss of precision by making the
 data type the same as the input::
 
     hdr = img.get_header()
     hdr.set_data_dtype(data.dtype)
-    img.write_data(data, fileobj)
+    img.to_filename(fname)
 
 '''
 
+import warnings
+
+
 class SpatialImage(object):
     _header_maker = dict
-    ''' Template class for lightweight image '''
+    ''' Template class for images '''
     def __init__(self, data, affine, header=None, extra=None):
         if extra is None:
             extra = {}
@@ -112,8 +134,17 @@
                 self.extra[key] = value
 
     @classmethod
-    def from_filespec(klass, filespec):
-        raise NotImplementedError
+    def from_filename(klass, filename):
+        files = klass.filespec_to_files(filename)
+        return klass.from_files(files)
+    
+    @classmethod
+    def from_filespec(klass, img, filespec):
+        warnings.warn('``from_filespec`` class method is deprecated\n'
+                      'Please use the ``from_filename`` class method '
+                      'instead',
+                      DeprecationWarning)
+        klass.from_filespec(filespec)
 
     def from_files(klass, files):
         raise NotImplementedError
@@ -125,8 +156,82 @@
     def filespec_to_files(filespec):
         raise NotImplementedError
     
-    def to_filespec(self, filespec):
-        raise NotImplementedError
+    def to_filename(self, filename):
+        ''' Write image to files implied by filename string
+
+        Paraameters
+        -----------
+        filename : str
+           filename to which to save image.  We will parse `filename`
+           with ``filespec_to_files`` to work out names for image,
+           header etc.
+
+        Returns
+        -------
+        None
+        '''
+        files = self.filespec_to_files(filename)
+        self.to_files(files)
+
+    def to_filespec(self, filename):
+        warnings.warn('``to_filespec`` is deprecated, please '
+                      'use ``to_filename`` instead',
+                      DeprecationWarning)
+        self.to_filename(filename)
 
     def to_files(self, files=None):
         raise NotImplementedError
+
+    @classmethod
+    def load(klass, filename):
+        return klass.from_filename(filename)
+
+    @classmethod
+    def save(klass, img, filename):
+        warnings.warn('``save`` class method is deprecated\n'
+                      'You probably want the ``to_filename`` instance '
+                      'method, or the module-level ``save`` function',
+                      DeprecationWarning)
+        klass.instance_to_filename(img, filename)
+
+    @classmethod
+    def instance_to_filename(klass, img, filename):
+        ''' Save `img` in our own format, to name implied by `filename`
+
+        This is a class method
+
+        Parameters
+        ----------
+        img : ``spatialimage`` instance
+           In fact, an object with the API of ``spatialimage`` -
+           specifically ``get_data``, ``get_affine``, ``get_header`` and
+           ``extra``.
+        filename : str
+           Filename, implying name to which to save image.
+        '''
+        img = klass.from_image(img)
+        img.to_filename(filename)
+        
+    @classmethod
+    def from_image(klass, img):
+        ''' Create new instance of own class from `img`
+
+        This is a class method
+        
+        Parameters
+        ----------
+        img : ``spatialimage`` instance
+           In fact, an object with the API of ``spatialimage`` -
+           specifically ``get_data``, ``get_affine``, ``get_header`` and
+           ``extra``.
+
+        Returns
+        -------
+        cimg : ``spatialimage`` instance
+           Image, of our own class
+        '''
+        return klass(img.get_data(),
+                     img.get_affine(),
+                     img.get_header(),
+                     img.extra)
+    

=== modified file 'nipy/io/imageformats/spm2analyze.py'
--- nipy/io/imageformats/spm2analyze.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/spm2analyze.py	2010-01-21 22:28:13 +0000
@@ -125,4 +125,4 @@
 
 
 load = Spm2AnalyzeImage.load
-save = Spm2AnalyzeImage.save
+save = Spm2AnalyzeImage.instance_to_filename

=== modified file 'nipy/io/imageformats/spm99analyze.py'
--- nipy/io/imageformats/spm99analyze.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/spm99analyze.py	2010-01-21 22:28:13 +0000
@@ -214,8 +214,8 @@
 class Spm99AnalyzeImage(analyze.AnalyzeImage):
     _header_maker = Spm99AnalyzeHeader
     @classmethod
-    def from_filespec(klass, filespec):
-        ret = super(Spm99AnalyzeImage, klass).from_filespec(filespec)
+    def from_filename(klass, filename):
+        ret = super(Spm99AnalyzeImage, klass).from_filename(filename)
         import scipy.io as sio
         matf = ret._files['mat']
         try:
@@ -274,4 +274,4 @@
 
 
 load = Spm99AnalyzeImage.load
-save = Spm99AnalyzeImage.save
+save = Spm99AnalyzeImage.instance_to_filename

=== modified file 'nipy/io/imageformats/testing/__init__.py'
--- nipy/io/imageformats/testing/__init__.py	2009-06-05 03:13:27 +0000
+++ nipy/io/imageformats/testing/__init__.py	2010-01-21 22:28:13 +0000
@@ -1,15 +1,18 @@
-''' Utilities to change context for nose.tools tests '''
-
-from os.path import join as pjoin, split as psplit, abspath
-
-import nose.tools as nt
-
-example_data_path = abspath(pjoin(psplit(__file__)[0], '..', 'tests', 'data'))
-
-assert_equal = lambda x, y: nt.assert_equal(x, y)
-assert_not_equal = lambda x, y: nt.assert_not_equal(x, y)
-assert_true = lambda x: nt.assert_true(x)
-assert_false = lambda x: nt.assert_false(x)
-
-def assert_raises(error, func, *args, **kwargs):
-    return nt.assert_raises(error, func, *args, **kwargs)
+''' Utilities for testing '''
+
+# Allow failed import of nose if not now running tests
+try:
+    import nose.tools as nt
+except ImportError:
+    pass
+else:
+    from lightunit import ParametricTestCase, parametric
+    # wrappers to change context for nose.tools tests '''
+    assert_equal = lambda x, y: nt.assert_equal(x, y)
+    assert_not_equal = lambda x, y: nt.assert_not_equal(x, y)
+    assert_true = lambda x: nt.assert_true(x)
+    assert_false = lambda x: nt.assert_false(x)
+
+    def assert_raises(error, func, *args, **kwargs):
+        return nt.assert_raises(error, func, *args, **kwargs)
+

=== added file 'nipy/io/imageformats/testing/_paramtestpy2.py'
--- nipy/io/imageformats/testing/_paramtestpy2.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/testing/_paramtestpy2.py	2010-01-21 22:28:13 +0000
@@ -0,0 +1,89 @@
+"""Implementation of the parametric test support for Python 2.x
+"""
+#-----------------------------------------------------------------------------
+# Imports
+#-----------------------------------------------------------------------------
+
+# Stdlib
+import unittest
+from compiler.consts import CO_GENERATOR
+
+#-----------------------------------------------------------------------------
+# Classes and functions
+#-----------------------------------------------------------------------------
+
+def isgenerator(func):
+    try:
+        return func.func_code.co_flags & CO_GENERATOR != 0
+    except AttributeError:
+        return False
+
+class ParametricTestCase(unittest.TestCase):
+    """Write parametric tests in normal unittest testcase form.
+
+    Limitations: the last iteration misses printing out a newline when running
+    in verbose mode.
+    """
+    def run_parametric(self, result, testMethod):
+        # But if we have a test generator, we iterate it ourselves
+        testgen = testMethod()
+        while True:
+            try:
+                # Initialize test
+                result.startTest(self)
+
+                # SetUp
+                try:
+                    self.setUp()
+                except KeyboardInterrupt:
+                    raise
+                except:
+                    result.addError(self, self._exc_info())
+                    return
+                # Test execution
+                ok = False
+                try:
+                    testgen.next()
+                    ok = True
+                except StopIteration:
+                    # We stop the loop
+                    break
+                except self.failureException:
+                    result.addFailure(self, self._exc_info())
+                except KeyboardInterrupt:
+                    raise
+                except:
+                    result.addError(self, self._exc_info())
+                # TearDown
+                try:
+                    self.tearDown()
+                except KeyboardInterrupt:
+                    raise
+                except:
+                    result.addError(self, self._exc_info())
+                    ok = False
+                if ok: result.addSuccess(self)
+
+            finally:
+                result.stopTest(self)
+
+    def run(self, result=None):
+        if result is None:
+            result = self.defaultTestResult()
+        testMethod = getattr(self, self._testMethodName)
+        # For normal tests, we just call the base class and return that
+        if isgenerator(testMethod):
+            return self.run_parametric(result, testMethod)
+        else:
+            return super(ParametricTestCase, self).run(result)
+
+
+def parametric(func):
+    """Decorator to make a simple function into a normal test via unittest."""
+
+    class Tester(ParametricTestCase):
+        test = staticmethod(func)
+
+    Tester.__name__ = func.__name__
+
+    return Tester

=== added file 'nipy/io/imageformats/testing/_paramtestpy3.py'
--- nipy/io/imageformats/testing/_paramtestpy3.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/testing/_paramtestpy3.py	2010-01-21 22:28:13 +0000
@@ -0,0 +1,62 @@
+"""Implementation of the parametric test support for Python 3.x.
+
+Thanks for the py3 version to Robert Collins, from the Testing in Python
+mailing list.
+"""
+#-----------------------------------------------------------------------------
+# Imports
+#-----------------------------------------------------------------------------
+
+# Stdlib
+import unittest
+from unittest import TestSuite
+
+#-----------------------------------------------------------------------------
+# Classes and functions
+#-----------------------------------------------------------------------------
+
+
+def isgenerator(func):
+    return hasattr(func,'_generator')
+
+
+class IterCallableSuite(TestSuite):
+   def __init__(self, iterator, adapter):
+       self._iter = iterator
+       self._adapter = adapter
+   def __iter__(self):
+       yield self._adapter(self._iter.__next__)
+
+class ParametricTestCase(unittest.TestCase):
+   """Write parametric tests in normal unittest testcase form.
+
+   Limitations: the last iteration misses printing out a newline when
+   running in verbose mode.
+   """
+
+   def run(self, result=None):
+       testMethod = getattr(self, self._testMethodName)
+       # For normal tests, we just call the base class and return that
+       if isgenerator(testMethod):
+           def adapter(next_test):
+               return unittest.FunctionTestCase(next_test,
+                                                self.setUp,
+                                                self.tearDown)
+
+           return IterCallableSuite(testMethod(),adapter).run(result)
+       else:
+           return super(ParametricTestCase, self).run(result)
+
+
+def parametric(func):
+   """Decorator to make a simple function into a normal test via
+unittest."""
+   # Hack, until I figure out how to write isgenerator() for python3!!
+   func._generator = True
+
+   class Tester(ParametricTestCase):
+       test = staticmethod(func)
+
+   Tester.__name__ = func.__name__
+
+   return Tester

=== added file 'nipy/io/imageformats/testing/lightunit.py'
--- nipy/io/imageformats/testing/lightunit.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/testing/lightunit.py	2010-01-21 22:28:13 +0000
@@ -0,0 +1,55 @@
+"""Lightweight testing that remains unittest-compatible.
+
+This module exposes decorators and classes to make lightweight testing possible
+in a manner similar to what nose allows, where standalone functions can be
+tests.  It also provides parametric test support that is vastly easier to use
+than nose's for debugging, because if a test fails, the stack under inspection
+is that of the test and not that of the test framework.
+
+- An @as_unittest decorator can be used to tag any normal parameter-less
+  function as a unittest TestCase.  Then, both nose and normal unittest will
+  recognize it as such.
+
+Authors
+-------
+
+- Fernando Perez <[e-mail]>
+"""
+
+#-----------------------------------------------------------------------------
+#  Copyright (C) 2009  The IPython Development Team
+#
+#  Distributed under the terms of the BSD License.  The full license is in
+#  the file COPYING, distributed as part of this software.
+#-----------------------------------------------------------------------------
+
+#-----------------------------------------------------------------------------
+# Imports
+#-----------------------------------------------------------------------------
+
+# Stdlib
+import sys
+import unittest
+
+# Our own
+import nosepatch
+
+if sys.version[0]=='2':
+    from _paramtestpy2 import ParametricTestCase, parametric
+else:
+    from _paramtestpy3 import ParametricTestCase, parametric
+
+#-----------------------------------------------------------------------------
+# Classes and functions
+#-----------------------------------------------------------------------------
+
+# Simple example of the basic idea
+def as_unittest(func):
+    """Decorator to make a simple function into a normal test via unittest."""
+    class Tester(unittest.TestCase):
+        def test(self):
+            func()
+
+    Tester.__name__ = func.__name__
+
+    return Tester

=== added file 'nipy/io/imageformats/testing/nosepatch.py'
--- nipy/io/imageformats/testing/nosepatch.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/testing/nosepatch.py	2010-01-21 22:28:13 +0000
@@ -0,0 +1,53 @@
+"""Monkeypatch nose to accept any callable as a method.
+
+By default, nose's ismethod() fails for static methods.
+Once this is fixed in upstream nose we can disable it.
+
+Note: merely importing this module causes the monkeypatch to be applied."""
+
+import unittest
+import nose.loader
+from inspect import ismethod, isfunction
+
+def getTestCaseNames(self, testCaseClass):
+    """Override to select with selector, unless
+    config.getTestCaseNamesCompat is True
+    """
+    if self.config.getTestCaseNamesCompat:
+        return unittest.TestLoader.getTestCaseNames(self, testCaseClass)
+
+    def wanted(attr, cls=testCaseClass, sel=self.selector):
+        item = getattr(cls, attr, None)
+        # MONKEYPATCH: replace this:
+        #if not ismethod(item):
+        #    return False
+        # return sel.wantMethod(item)
+        # With:
+        if ismethod(item):
+            return sel.wantMethod(item)
+        # static method or something. If this is a static method, we
+        # can't get the class information, and we have to treat it
+        # as a function.  Thus, we will miss things like class
+        # attributes for test selection
+        if isfunction(item):
+            return sel.wantFunction(item)
+        return False
+        # END MONKEYPATCH
+    
+    cases = filter(wanted, dir(testCaseClass))
+    for base in testCaseClass.__bases__:
+        for case in self.getTestCaseNames(base):
+            if case not in cases:
+                cases.append(case)
+    # add runTest if nothing else picked
+    if not cases and hasattr(testCaseClass, 'runTest'):
+        cases = ['runTest']
+    if self.sortTestMethodsUsing:
+        cases.sort(self.sortTestMethodsUsing)
+    return cases
+
+
+##########################################################################
+# Apply monkeypatch here
+nose.loader.TestLoader.getTestCaseNames = getTestCaseNames
+##########################################################################

=== modified file 'nipy/io/imageformats/tests/test_funcs.py'
--- nipy/io/imageformats/tests/test_funcs.py	2009-12-12 19:52:25 +0000
+++ nipy/io/imageformats/tests/test_funcs.py	2010-01-21 22:28:13 +0000
@@ -6,10 +6,15 @@
 
 import nipy.io.imageformats as nf
 
-from nipy.io.imageformats.funcs import concat_images
+from nipy.io.imageformats.funcs import concat_images, as_closest_canonical, OrientationError
+from nipy.io.imageformats.nifti1 import Nifti1Image
+
+from nipy.testing import parametric
 
 from numpy.testing import assert_array_equal
-from nose.tools import assert_true, assert_equal, assert_raises
+from nose.tools import (assert_true, assert_false,
+                        assert_equal, assert_raises)
+
 
 def test_concat():
     shape = (1,2,5)
@@ -28,5 +33,28 @@
     img2 = nf.Nifti1Image(data1.T, affine)
     yield assert_raises, ValueError, concat_images, [img0, img2]
 
-    
-    
+
+@parametric
+def test_closest_canonical():
+    arr = np.arange(24).reshape((2,3,4,1))
+    # no funky stuff, returns same thing
+    img = Nifti1Image(arr, np.eye(4))
+    xyz_img = as_closest_canonical(img)
+    yield assert_true(img is xyz_img)
+    # a axis flip
+    img = Nifti1Image(arr, np.diag([-1,1,1,1]))
+    xyz_img = as_closest_canonical(img)
+    yield assert_false(img is xyz_img)
+    out_arr = xyz_img.get_data()
+    yield assert_array_equal(out_arr, np.flipud(arr))
+    # no error for enforce_diag in this case
+    xyz_img = as_closest_canonical(img, True)
+    # but there is if the affine is not diagonal
+    aff = np.eye(4)
+    aff[0,1] = 0.1
+    # although it's more or less canonical already
+    img = Nifti1Image(arr, aff)
+    xyz_img = as_closest_canonical(img)
+    yield assert_true(img is xyz_img)
+    # it's still not diagnonal
+    yield assert_raises(OrientationError, as_closest_canonical, img, True)

=== modified file 'nipy/io/imageformats/tests/test_nifti1.py'
--- nipy/io/imageformats/tests/test_nifti1.py	2009-09-11 01:35:49 +0000
+++ nipy/io/imageformats/tests/test_nifti1.py	2010-01-21 22:28:13 +0000
@@ -270,7 +270,7 @@
     for ext in ('.gz', '.bz2'):
         try:
             _, fname = tempfile.mkstemp('.nii' + ext)
-            img.to_filespec(fname)
+            img.to_filename(fname)
             img3 = Nifti1Image.load(fname)
             yield assert_true, isinstance(img3, img.__class__)
             yield assert_array_equal, img3.get_data(), data

=== added file 'nipy/io/imageformats/tests/test_orientations.py'
--- nipy/io/imageformats/tests/test_orientations.py	1970-01-01 00:00:00 +0000
+++ nipy/io/imageformats/tests/test_orientations.py	2010-01-21 22:28:13 +0000
@@ -0,0 +1,219 @@
+''' Testing for orientations module '''
+
+import numpy as np
+
+from nose.tools import assert_true, assert_equal, assert_raises
+
+from numpy.testing import assert_array_equal, assert_array_almost_equal
+
+from nipy.io.imageformats.orientations import (io_orientation, orientation_affine,
+                                  flip_axis, _ornt_to_affine,
+                                  apply_orientation, OrientationError)
+
+from nipy.io.imageformats.testing import parametric
+
+
+IN_ARRS = [np.eye(4),
+           [[0,0,1,0],
+            [0,1,0,0],
+            [1,0,0,0],
+            [0,0,0,1]],
+           [[0,1,0,0],
+            [0,0,1,0],
+            [1,0,0,0],
+            [0,0,0,1]],
+           [[3,1,0,0],
+            [1,3,0,0],
+            [0,0,1,0],
+            [0,0,0,1]],
+           [[1,3,0,0],
+            [3,1,0,0],
+            [0,0,1,0],
+            [0,0,0,1]],
+           ]
+
+OUT_ORNTS = [[[0,1],
+              [1,1],
+              [2,1]],
+             [[2,1],
+              [1,1],
+              [0,1]],
+             [[2,1],
+              [0,1],
+              [1,1]],
+             [[0,1],
+              [1,1],
+              [2,1]],
+             [[1,1],
+              [0,1],
+              [2,1]],
+            ]
+
+IN_ARRS = IN_ARRS + [[[np.cos(np.pi/6+i*np.pi/2),np.sin(np.pi/6+i*np.pi/2),0,0], 
+                      [-np.sin(np.pi/6+i*np.pi/2),np.cos(np.pi/6+i*np.pi/2),0,0],
+                      [0,0,1,0],
+                      [0,0,0,1]] for i in range(4)]
+
+OUT_ORNTS = OUT_ORNTS + [[[0,1],
+                          [1,1],
+                          [2,1]],
+                         [[1,-1],
+                          [0,1],
+                          [2,1]],
+                         [[0,-1],
+                          [1,-1],
+                          [2,1]],
+                         [[1,1],
+                          [0,-1],
+                          [2,1]]
+                         ]
+
+                         
+IN_ARRS = [np.array(arr) for arr in IN_ARRS]
+OUT_ORNTS = [np.array(ornt) for ornt in OUT_ORNTS]
+
+
+def same_transform(taff, ornt, shape):
+    # Applying transformations implied by `ornt` to a made-up array
+    # ``arr`` of shape `shape`, results in ``t_arr``. When the point
+    # indices from ``arr`` are transformed by (the inverse of) `taff`,
+    # and we index into ``t_arr`` with these transformed points, then we
+    # should get the same values as we would from indexing into arr with
+    # the untransformed points. 
+    shape = np.array(shape)
+    size = np.prod(shape)
+    arr = np.arange(size).reshape(shape)
+    # apply ornt transformations
+    t_arr = apply_orientation(arr, ornt)
+    # get all point indices in arr
+    i,j,k = shape
+    arr_pts = np.mgrid[:i,:j,:k].reshape((3,-1))
+    # inverse of taff takes us from point index in arr to point index in
+    # t_arr
+    itaff = np.linalg.inv(taff)
+    # apply itaff so that points indexed in t_arr should correspond 
+    o2t_pts = np.dot(itaff[:3,:3], arr_pts) + itaff[:3,3][:,None]
+    assert np.allclose(np.round(o2t_pts), o2t_pts)
+    # fancy index out the t_arr values
+    vals = t_arr[list(o2t_pts.astype('i'))]
+    return np.all(vals == arr.ravel())
+
+
+@parametric
+def test_apply():
+    # most tests are in ``same_transform`` above, via the
+    # test_io_orientations. 
+    a = np.arange(24).reshape((2,3,4))
+    # Test 4D
+    t_arr = apply_orientation(a[:,:,:,None], ornt)
+    yield assert_equal(t_arr.ndim, 4)
+    # Orientation errors
+    yield assert_raises(OrientationError,
+                        apply_orientation,
+                        a[:,:,1], ornt)
+    yield assert_raises(OrientationError,
+                        apply_orientation,
+                        a,
+                        [[0,1],[np.nan,np.nan],[2,1]])
+                        
+
+@parametric
+def test_flip_axis():
+    a = np.arange(24).reshape((2,3,4))
+    yield assert_array_equal(
+        flip_axis(a),
+        np.flipud(a))
+    yield assert_array_equal(
+        flip_axis(a, axis=0),
+        np.flipud(a))
+    yield assert_array_equal(
+        flip_axis(a, axis=1),
+        np.fliplr(a))
+    # check accepts array-like
+    yield assert_array_equal(
+        flip_axis(a.tolist(), axis=0),
+        np.flipud(a))
+    # third dimension
+    b = a.transpose()
+    b = np.flipud(b)
+    b = b.transpose()
+    yield assert_array_equal(flip_axis(a, axis=2), b)
+
+
+@parametric
+def test_io_orientation():
+    shape = (2,3,4)
+    for in_arr, out_ornt in zip(IN_ARRS, OUT_ORNTS):
+        ornt = io_orientation(in_arr)
+        yield assert_array_equal(ornt, out_ornt)
+        taff = orientation_affine(ornt, shape)
+        yield assert_true(same_transform(taff, ornt, shape))
+        for axno in range(3):
+            arr = in_arr.copy()
+            ex_ornt = out_ornt.copy()
+            # flip the input axis in affine
+            arr[:,axno] *= -1
+            # check that result shows flip
+            ex_ornt[axno, 1] *= -1
+            ornt = io_orientation(arr)
+            yield assert_array_equal(ornt, ex_ornt)
+            taff = orientation_affine(ornt, shape)
+            yield assert_true(same_transform(taff, ornt, shape))
+           
+
+def test_drop_coord():
+    # given a 5x4 affine from slicing an fmri,
+    # the orientations code should easily reorder and drop the t
+    # axis
+
+    # this affine has output coordinates '-y','z','x' and is at t=16
+    sliced_fmri_affine = np.array([[0,-1,0,3],
+                                   [0,0,2,5],
+                                   [3,0,0,4],
+                                   [0,0,0,16],
+                                   [0,0,0,1]])
+    ornt = io_orientation(sliced_fmri_affine)
+    affine_that_drops_t_reorders_and_flips = _ornt_to_affine(ornt)
+    final_affine = np.dot(affine_that_drops_t_reorders_and_flips, 
+                          sliced_fmri_affine)
+    # the output will be diagonal
+    # with the 'y' row having been flipped and the 't' row dropped
+    assert_array_equal(final_affine, 
+                       np.array([[3,0,0,4],
+                                 [0,1,0,-3],
+                                 [0,0,2,5],
+                                 [0,0,0,1]]))
+
+
+@parametric
+def test_ornt_to_affine():
+    # this orientation indicates that the first output
+    # axis of the affine is closest to the vector [0,0,-1],
+    # the last is closest to [1,0,0] and 
+    # that the y coordinate ([0,1,0]) is dropped
+    ornt = [[2,-1],
+            [np.nan,np.nan],
+            [0,1]]
+    # the reordering/flipping is represented by an affine that 
+    # takes the 3rd output coordinate and maps it to the
+    # first, takes the 3rd, maps it to first and flips it
+    A = np.array([[0,0,-1,0],
+                  [1,0,0,0],
+                  [0,0,0,1]])
+    yield assert_array_equal(A, _ornt_to_affine(ornt))
+    # a more complicated example. only the 1st, 3rd and 6th
+    # coordinates appear in the output
+    ornt = [[3,-1],
+            [np.nan,np.nan],
+            [0,1],
+            [np.nan,np.nan],
+            [np.nan,np.nan],
+            [1,-1]]
+    B = np.array([[0,0,0,-1,0,0,0],
+                  [1,0,0,0,0,0,0],
+                  [0,-1,0,0,0,0,0],
+                  [0,0,0,0,0,0,1]])
+    yield assert_array_equal(B, _ornt_to_affine(ornt))
+
+
+    

=== modified file 'nipy/io/nifti_ref.py'
--- nipy/io/nifti_ref.py	2009-11-25 17:08:03 +0000
+++ nipy/io/nifti_ref.py	2010-01-21 22:28:12 +0000
@@ -49,7 +49,7 @@
 
 FIXME: (finish this thought.... Are we going to open NIFTI files with
 NIFTI input coordinates?)  For this reason, we have to consider whether
-we should transpose the memmap from pynifti.
+we should transpose the memmap from nifti input.
 """
 import warnings
 

=== modified file 'nipy/modalities/fmri/fmristat/tests/test_FIAC.py'
--- nipy/modalities/fmri/fmristat/tests/test_FIAC.py	2009-12-03 21:12:30 +0000
+++ nipy/modalities/fmri/fmristat/tests/test_FIAC.py	2010-01-21 22:28:13 +0000
@@ -306,10 +306,9 @@
     i = interp1d(t, formula.vectorize(col)(t))
     return formula.aliased_function(name, i)(formula.t)
 
+
 def test_agreement():
-    """
-    The test: does Protocol manage to recreate the design of fMRIstat?
-    """
+    # The test: does Protocol manage to recreate the design of fMRIstat?
     for design_type in ['event', 'block']:
         dd = D[design_type]
 
@@ -319,20 +318,18 @@
                 yield niptest.assert_true, np.greater(cmax, 0.999)
 
 
-
 def test_event_design():
-
-	block = csv2rec(StringIO(altdescr['block']))
-	event = csv2rec(StringIO(altdescr['event']))
-	t = np.arange(191)*2.5+1.25
-			
-	bkeep = np.not_equal((np.arange(block.time.shape[0])) % 6, 0)
-	ekeep = np.greater(np.arange(event.time.shape[0]), 0)
-
-	# Even though there is a FIAC block experiment
-	# the design is represented as an event design
-	# with the same event repeated several times in a row...
-
-	Xblock, cblock = design.event_design(block[bkeep], t, hrfs=delay.spectral)
-	Xevent, cevent = design.event_design(event[ekeep], t, hrfs=delay.spectral)
-
+    block = csv2rec(StringIO(altdescr['block']))
+    event = csv2rec(StringIO(altdescr['event']))
+    t = np.arange(191)*2.5+1.25
+    
+    bkeep = np.not_equal((np.arange(block.time.shape[0])) % 6, 0)
+    ekeep = np.greater(np.arange(event.time.shape[0]), 0)
+    
+    # Even though there is a FIAC block experiment
+    # the design is represented as an event design
+    # with the same event repeated several times in a row...
+    
+    Xblock, cblock = design.event_design(block[bkeep], t, hrfs=delay.spectral)
+    Xevent, cevent = design.event_design(event[ekeep], t, hrfs=delay.spectral)
+    

=== modified file 'nipy/modalities/fmri/pca.py'
--- nipy/modalities/fmri/pca.py	2009-12-16 18:35:13 +0000
+++ nipy/modalities/fmri/pca.py	2010-01-21 22:28:13 +0000
@@ -13,7 +13,7 @@
 """
 
 import numpy as np
-import numpy.linalg as npl
+import scipy.linalg as spl
 from nipy.fixes.scipy.stats.models.utils import pos_recipr
 
 
@@ -25,36 +25,39 @@
     ----------
     data : ndarray-like (np.float)
        The array on which to perform PCA over axis `axis` (below)
-    axis : int
+    axis : int, optional
        The axis over which to perform PCA (axis identifying
        observations).  Default is 0 (first)
-    mask : ndarray-like (np.bool)
+    mask : ndarray-like (np.bool), optional
        An optional mask, should have shape given by data axes, with
-       `axis` removed, i.e.: ``s = data.shape; s.pop(axis); msk_shape=s``
-    ncomp : {None, int}
+       `axis` removed, i.e.: ``s = data.shape; s.pop(axis);
+       msk_shape=s``
+    ncomp : {None, int}, optional
        How many component basis projections to return. If ncomp is None
        (the default) then the number of components is given by the
        calculated rank of the data, after applying `design_keep`,
        `design_resid` and `tol_ratio` below.  We always return all the
        basis vectors and percent variance for each component; `ncomp`
        refers only to the number of basis_projections returned.
-    standardize : bool
-       Standardize so each time series has same error-sum-of-squares?
-    design_keep : None or ndarray
+    standardize : bool, optional
+       If True, standardize so each time series (after application of
+       `design_keep` and `design_resid`) has the same standard
+       deviation, as calculated by the ``np.std`` function.
+    design_keep : None or ndarray, optional
        Data is projected onto the column span of design_keep.
        None (default) equivalent to ``np.identity(data.shape[axis])``
-    design_resid : str or None or ndarray
+    design_resid : str or None or ndarray, optional
        After projecting onto the column span of design_keep, data is
        projected perpendicular to the column span of this matrix.  If
        None, we do no such second projection.  If a string 'mean', then
        the mean of the data is removed, equivalent to passing a column
        vector matrix of 1s.
-    tol_ratio : float
+    tol_ratio : float, optional
        If ``XZ`` is the vector of singular values of the projection
        matrix from `design_keep` and `design_resid`, and S are the
        singular values of ``XZ``, then `tol_ratio` is the value used to
-       calculate the effective rank of the projection, as in ``rank = ((S
-       / S.max) > tol_ratio).sum()``
+       calculate the effective rank of the projection of the design, as
+       in ``rank = ((S / S.max) > tol_ratio).sum()``
 
     Returns
     -------
@@ -65,23 +68,24 @@
     
        `results` has keys:
 
-       * ``basis_vectors``: time series, shape (data.shape[axis], L)
+       * ``basis_vectors``: series over `axis`, shape (data.shape[axis], L) -
+          the eigenvectors of the PCA
        * ``pcnt_var``: percent variance explained by component, shape
           (L,)
-       * ``basis_projections``: PCA components, with components varying over axis
-          `axis`; thus shape given by: ``s = list(data.shape); s[axis] =
-          ncomp``
+       * ``basis_projections``: PCA components, with components varying
+          over axis `axis`; thus shape given by: ``s = list(data.shape);
+          s[axis] = ncomp``
        * ``axis``: axis over which PCA has been performed.
 
-     Notes
-     -----
-     See ``pca_image.m`` from fmristat for Keith Worsley's code on which
-     some of this is based.
+    Notes
+    -----
+    See ``pca_image.m`` from ``fmristat`` for Keith Worsley's code on
+    which some of this is based.
 
-     See: http://en.wikipedia.org/wiki/Principal_component_analysis for
-     some inspiration for naming - particularly 'basis_vectors' and
-     'basis_projections'
-     """
+    See: http://en.wikipedia.org/wiki/Principal_component_analysis for
+    some inspiration for naming - particularly 'basis_vectors' and
+    'basis_projections'
+    """
     data = np.asarray(data)
     # We roll the PCA axis to be first, for convenience
     if axis is None:
@@ -89,7 +93,6 @@
     data = np.rollaxis(data, axis)
     if mask is not None:
         mask = np.asarray(mask)
-    nbasis_projections = data.shape[0]
     if design_resid == 'mean':
         # equivalent to: design_resid = np.ones((data.shape[0], 1))
         def project_resid(Y):
@@ -97,9 +100,20 @@
     elif design_resid is None:
         def project_resid(Y): return Y
     else: # matrix passed, we hope
-        pinv_design_resid = npl.pinv(design_resid)
+        projector = np.dot(design_resid, spl.pinv(design_resid))
         def project_resid(Y):
-            return Y - np.dot(np.dot(design_resid, pinv_design_resid), Y)
+            return Y - np.dot(projector, Y)
+    if standardize:
+        def standardize_from(arr, std_source):
+            # modifies array in place
+            resid = project_resid(std_source)
+            rstd = np.sqrt(np.square(resid).sum(axis=0) / resid.shape[0])
+            # positive 1/rstd
+            rstd_half = np.where(rstd<=0, 0, 1. / rstd)
+            arr *= rstd_half
+            return arr
+    else:
+        standardize_from = None
     """
     Perform the computations needed for the PCA.  This stores the
     covariance/correlation matrix of the data in the attribute 'C'.  The
@@ -111,20 +125,20 @@
     to column space of design_resid.
     """
     if design_keep is None:
-        X = np.eye(nbasis_projections)
+        X = np.eye(data.shape[0])
     else:
-        X = np.dot(design_keep, npl.pinv(design_keep))
+        X = np.dot(design_keep, spl.pinv(design_keep))
     XZ = project_resid(X)
-    UX, SX, VX = npl.svd(XZ, full_matrices=0)
+    UX, SX, VX = spl.svd(XZ, full_matrices=0)
     # The matrix UX has orthonormal columns and represents the
     # final "column space" that the data will be projected onto.
     rank = (SX/SX.max() > tol_ratio).sum()
     UX = UX[:,range(rank)].T
     # calculate covariance matrix
-    C  = _get_covariance(data, UX, standardize, project_resid, mask)
+    C  = _get_covariance(data, UX, standardize_from, mask)
     # find the eigenvalues D and eigenvectors Vs of the covariance
     # matrix
-    D, Vs = npl.eigh(C)
+    D, Vs = spl.eigh(C)
     # sort both in descending order of eigenvalues
     order = np.argsort(-D)
     D = D[order]
@@ -136,7 +150,7 @@
     if ncomp is None:
         ncomp = rank
     subVX = basis_vectors[:ncomp]
-    out = _get_basis_projections(data, subVX)
+    out = _get_basis_projections(data, subVX, standardize_from)
     # Roll PCA image axis back to original position in data array
     if axis < 0:
         axis += data.ndim
@@ -147,29 +161,33 @@
             'axis': axis}
 
 
-def _get_covariance(data, UX, standardize, project_resid, mask):
+def _get_covariance(data, UX, standardize_from, mask):
+    # number of points in PCA dimension
     rank = UX.shape[0]
+    n_pts = data.shape[0]
     C = np.zeros((rank, rank))
+    # loop over next dimension to save memory
     for i in range(data.shape[1]):
-        Y = data[:,i].reshape((data.shape[0], -1))
+        Y = data[:,i].reshape((n_pts, -1))
         # project data into required space
         YX = np.dot(UX, Y)
-        if standardize:
-            S2 = (project_resid(Y)**2).sum(0)
-            Smhalf = pos_recipr(np.sqrt(S2)); del(S2)
-            YX *= Smhalf
+        if standardize_from is not None:
+            YX = standardize_from(YX, Y)
         if mask is not None:
+            # weight data with mask.  Usually the weights will be 0,1
             YX = YX * np.nan_to_num(mask[i].reshape(Y.shape[1]))
         C += np.dot(YX, YX.T)
     return C
 
 
-def _get_basis_projections(data, subVX):
+def _get_basis_projections(data, subVX, standardize_from):
     ncomp = subVX.shape[0]
     out = np.empty((ncomp,) + data.shape[1:], np.float)
     for i in range(data.shape[1]):
         Y = data[:,i].reshape((data.shape[0], -1))
         U = np.dot(subVX, Y)
+        if standardize_from is not None:
+           U = standardize_from(U, Y)
         U.shape = (U.shape[0],) + data.shape[2:]
         out[:,i] = U
     return out

=== modified file 'nipy/modalities/fmri/tests/test_pca.py'
--- nipy/modalities/fmri/tests/test_pca.py	2009-12-14 22:01:46 +0000
+++ nipy/modalities/fmri/tests/test_pca.py	2010-01-21 22:28:13 +0000
@@ -36,13 +36,22 @@
     return recond
 
 
+def root_mse(arr, axis=0):
+    return np.sqrt(np.square(arr).sum(axis=axis) / arr.shape[axis])
+
+
+@parametric
 def test_input_effects():
     ntotal = data['nimages'] - 1
     # return full rank - mean PCA over last axis
     p = pca(data['fmridata'], -1)
-    yield assert_equal, p['basis_vectors'].shape, (data['nimages'], ntotal)
-    yield assert_equal, p['basis_projections'].shape, data['mask'].shape + (ntotal,)
-    yield assert_equal, p['pcnt_var'].shape, (ntotal,)
+    yield assert_equal(
+        p['basis_vectors'].shape,
+        (data['nimages'], ntotal))
+    yield assert_equal(
+        p['basis_projections'].shape,
+        data['mask'].shape + (ntotal,))
+    yield assert_equal(p['pcnt_var'].shape, (ntotal,))
     # Reconstructed data lacks only mean
     rarr = reconstruct(p['basis_vectors'], p['basis_projections'], -1)
     rarr = rarr + data['fmridata'].mean(-1)[...,None]
@@ -51,19 +60,19 @@
     arr = np.rollaxis(arr, -1)
     pr = pca(arr)
     out_arr = np.rollaxis(pr['basis_projections'], 0, 4)
-    yield assert_array_equal, out_arr, p['basis_projections']
-    yield assert_array_equal, p['basis_vectors'], pr['basis_vectors']
-    yield assert_array_equal, p['pcnt_var'], pr['pcnt_var']
+    yield assert_array_equal(out_arr, p['basis_projections'])
+    yield assert_array_equal(p['basis_vectors'], pr['basis_vectors'])
+    yield assert_array_equal(p['pcnt_var'], pr['pcnt_var'])
     # Check axis None raises error
-    yield assert_raises, ValueError, pca, data['fmridata'], None
+    yield assert_raises(ValueError, pca, data['fmridata'], None)
 
 
 @parametric
 def test_diagonality():
-    # basis_projections are not diagonal, unless not standarized
-    p = pca(data['fmridata'], -1)
-    yield assert_false(diagonal_covariance(p['basis_projections'], -1))
-    pns = pca(data['fmridata'], -1, standardize=False)
+    # basis_projections are diagonal, whether standarized or not
+    p = pca(data['fmridata'], -1) # standardized
+    yield assert_true(diagonal_covariance(p['basis_projections'], -1))
+    pns = pca(data['fmridata'], -1, standardize=False) # not 
     yield assert_true(diagonal_covariance(pns['basis_projections'], -1))
 
 
@@ -74,6 +83,7 @@
     return np.allclose(aTa, np.diag(np.diag(aTa)), atol=1e-6)
 
 
+@parametric
 def test_2D():
     # check that a standard 2D PCA works too
     M = 100
@@ -83,26 +93,35 @@
     p = pca(data)
     ts = p['basis_vectors']
     imgs = p['basis_projections']
-    yield assert_equal, ts.shape, (M, L)
-    yield assert_equal, imgs.shape, (L, N)
-    rimgs = reconstruct(ts, imgs) + data.mean(0)[None,...]
-    yield assert_array_almost_equal, rimgs, data
-    # if standardize is set, projections will not be covariance diagonal
-    yield assert_false, diagonal_covariance(imgs)
-    # but will if not
+    yield assert_equal(ts.shape, (M, L))
+    yield assert_equal(imgs.shape, (L, N))
+    rimgs = reconstruct(ts, imgs)
+    # add back the sqrt MSE, because we standardized
+    data_mean = data.mean(0)[None,...]
+    demeaned = data - data_mean
+    rmse = root_mse(demeaned, axis=0)[None,...]
+    # also add back the mean
+    yield assert_array_almost_equal((rimgs * rmse) + data_mean, data)
+    # if standardize is set, or not, covariance is diagonal
+    yield assert_true(diagonal_covariance(imgs))
     p = pca(data, standardize=False)
     imgs = p['basis_projections']
-    yield assert_true, diagonal_covariance(imgs)
+    yield assert_true(diagonal_covariance(imgs))
     
 
+@parametric
 def test_PCAMask():
     ntotal = data['nimages'] - 1
     ncomp = 5
     p = pca(data['fmridata'], -1, data['mask'], ncomp=ncomp)
-    yield assert_equal, p['basis_vectors'].shape, (data['nimages'], ntotal)
-    yield assert_equal, p['basis_projections'].shape, data['mask'].shape + (ncomp,)
-    yield assert_equal, p['pcnt_var'].shape, (ntotal,)
-    yield assert_almost_equal, p['pcnt_var'].sum(), 100.
+    yield assert_equal(
+        p['basis_vectors'].shape,
+        (data['nimages'], ntotal))
+    yield assert_equal(
+        p['basis_projections'].shape,
+        data['mask'].shape + (ncomp,))
+    yield assert_equal(p['pcnt_var'].shape, (ntotal,))
+    yield assert_almost_equal(p['pcnt_var'].sum(), 100.)
 
 
 def test_PCAMask_nostandardize():
@@ -151,6 +170,7 @@
     yield assert_almost_equal, p['pcnt_var'].sum(), 100.
 
 
+@parametric
 def test_resid():
     # Data is projected onto k=10 dimensional subspace then has its mean
     # removed.  Should still have rank 10.
@@ -159,17 +179,23 @@
     ntotal = k
     X = np.random.standard_normal((data['nimages'], k))
     p = pca(data['fmridata'], -1, ncomp=ncomp, design_resid=X)
-    yield assert_equal, p['basis_vectors'].shape, (data['nimages'], ntotal)
-    yield assert_equal, p['basis_projections'].shape, data['mask'].shape + (ncomp,)
-    yield assert_equal, p['pcnt_var'].shape, (ntotal,)
-    yield assert_almost_equal, p['pcnt_var'].sum(), 100.
+    yield assert_equal(
+        p['basis_vectors'].shape,
+        (data['nimages'], ntotal))
+    yield assert_equal(
+        p['basis_projections'].shape,
+        data['mask'].shape + (ncomp,))
+    yield assert_equal(p['pcnt_var'].shape, (ntotal,))
+    yield assert_almost_equal(p['pcnt_var'].sum(), 100.)
     # if design_resid is None, we do not remove the mean, and we get
     # full rank from our data
     p = pca(data['fmridata'], -1, design_resid=None)
     rank = p['basis_vectors'].shape[1]
-    yield assert_equal, rank, data['nimages']
+    yield assert_equal(rank, data['nimages'])
     rarr = reconstruct(p['basis_vectors'], p['basis_projections'], -1)
-    yield assert_array_almost_equal, rarr, data['fmridata']
+    # add back the sqrt MSE, because we standardized
+    rmse = root_mse(data['fmridata'], axis=-1)[...,None]
+    yield assert_array_almost_equal(rarr * rmse, data['fmridata'])
 
 
 def test_both():

=== modified file 'nipy/neurospin/register/benchmarks/bench_register.py'
--- nipy/neurospin/register/benchmarks/bench_register.py	2009-07-11 17:35:34 +0000
+++ nipy/neurospin/register/benchmarks/bench_register.py	2010-01-21 22:28:13 +0000
@@ -39,7 +39,7 @@
     dt2 = time.clock()-t0
     assert_almost_equal(I1, I2)
     print('3d array resampling')
-    print('  using neuroimaging.neurospin: %f sec' % dt1)
+    print('  using nipy.neurospin: %f sec' % dt1)
     print('  using scipy.ndimage: %f sec' % dt2)
 
 

=== modified file 'setup.py'
--- setup.py	2009-12-17 23:36:20 +0000
+++ setup.py	2010-01-21 22:28:12 +0000
@@ -18,11 +18,8 @@
     #    'nipy.core.image')
     # Robert Kern recommends setting quiet=True on the numpy list, stating
     # these messages are probably only used in debugging numpy distutils.
-
     config.get_version('nipy/version.py') # sets config.version
-
     config.add_subpackage('nipy', 'nipy')
-
     return config
 
 ################################################################################
@@ -41,28 +38,76 @@
 
 
 # Dependency checks
-def package_check(pkg_name, version=None, checker=LooseVersion):
-    msg = 'Missing package: %s, you might get run-time errors' % pkg_name 
+def package_check(pkg_name, version=None,
+                  optional=False,
+                  checker=LooseVersion,
+                  version_getter=None,
+                  ):
+    ''' Check if package `pkg_name` is present, and correct version
+
+    Parameters
+    ----------
+    pkg_name : str
+       name of package as imported into python
+    version : {None, str}, optional
+       minimum version of the package that we require. If None, we don't
+       check the version.  Default is None
+    optional : {False, True}, optional
+       If False, raise error for absent package or wrong version;
+       otherwise warn
+    checker : callable, optional
+       callable with which to return comparable thing from version
+       string.  Default is ``distutils.version.LooseVersion``
+    version_getter : {None, callable}:
+       Callable that takes `pkg_name` as argument, and returns the
+       package version string - as in::
+       
+          ``version = version_getter(pkg_name)``
+
+       If None, equivalent to::
+
+          mod = __import__(pkg_name); version = mod.__version__``
+    '''
+    if version_getter is None:
+        def version_getter(pkg_name):
+            mod = __import__(pkg_name)
+            return mod.__version__
     try:
         mod = __import__(pkg_name)
     except ImportError:
-        log.warn(msg)
-        return 
+        if not optional:
+            raise RuntimeError('Cannot import package "%s" '
+                               '- is it installed?' % pkg_name)
+        log.warn('Missing optional package "%s"; '
+                 'you may get run-time errors' % pkg_name)
+        return
     if not version:
         return
-    msg += ' >= %s' % version
     try:
-        have_version = mod.__version__
+        have_version = version_getter(pkg_name)
     except AttributeError:
         raise RuntimeError('Cannot find version for %s' % pkg_name)
     if checker(have_version) < checker(version):
-        raise RuntimeError(msg)
-
-
-# Soft dependency checking
-package_check('sympy', '0.6.4')
+        v_msg = 'You have version %s of package "%s"' \
+            ' but we need version >= %s' % (
+            have_version,
+            pkg_name,
+            version,
+            )
+        if optional:
+            log.warn(v_msg + '; you may get run-time errors')
+        else:
+            raise RuntimeError(v_msg)
+
+
+# Hard and soft dependency checking
 package_check('scipy', '0.5')
-
+package_check('sympy', '0.6.6')
+def _mayavi_version(pkg_name):
+    from enthought.mayavi import version
+    return version.version
+package_check('mayavi', '3.0', optional=True,
+              version_getter=_mayavi_version)
     
 ################################################################################
 # Import the documentation building classes.