def cortex_to_image_interpolation(obj, mask=None, affine=None, method='<STR_LIT>', shape=None):
<EOL>if mask is None:<EOL><INDENT>if shape is None and affine is None:<EOL><INDENT>shape = (<NUM_LIT>,<NUM_LIT>,<NUM_LIT>)<EOL>affine = ([[-<NUM_LIT:1>,<NUM_LIT:0>,<NUM_LIT:0>],[<NUM_LIT:0>,<NUM_LIT:0>,-<NUM_LIT:1>],[<NUM_LIT:0>,<NUM_LIT:1>,<NUM_LIT:0>]], [<NUM_LIT>,<NUM_LIT>,<NUM_LIT>])<EOL><DEDENT>elif shape is None:<EOL><INDENT>shape = (<NUM_LIT>,<NUM_LIT>,<NUM_LIT>)<EOL><DEDENT>elif affine is None:<EOL><INDENT>affine = ([[-<NUM_LIT:1>,<NUM_LIT:0>,<NUM_LIT:0>],[<NUM_LIT:0>,<NUM_LIT:0>,-<NUM_LIT:1>],[<NUM_LIT:0>,<NUM_LIT:1>,<NUM_LIT:0>]], np.asarray(shape) / <NUM_LIT:2>)<EOL><DEDENT>mask = obj.gray_indices if isinstance(obj, Subject) else '<STR_LIT:all>'<EOL><DEDENT>if pimms.is_matrix(mask):<EOL><INDENT>if not is_tuple(mask) or len(mask) != <NUM_LIT:3>:<EOL><INDENT>mask = np.asarray(mask, dtype=np.int)<EOL>mask = tuple(mask.T if mask.shape[<NUM_LIT:0>] != <NUM_LIT:3> else mask)<EOL><DEDENT><DEDENT>elif isinstance(mask, colls.Set):<EOL><INDENT>mask = np.asarray(list(mask)).T<EOL>tmp = np.full(shape, False)<EOL>tmp[tuple(mask)] = True<EOL>mask = np.where(tmp)<EOL><DEDENT>elif pimms.is_array(mask, dims=<NUM_LIT:3>):<EOL><INDENT>if shape is None: shape = mask.shape<EOL>mask = np.where(mask)<EOL><DEDENT>elif pimms.is_str(mask) and mask.lower() == '<STR_LIT:all>':<EOL><INDENT>if shape is None: shape = (<NUM_LIT>,<NUM_LIT>,<NUM_LIT>)<EOL>mask = np.where(np.ones(shape))<EOL><DEDENT>else:<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>if shape is None:  shape  = np.asarray(np.ceil(np.max(mask, axis=<NUM_LIT:0>)), dtype=np.int)<EOL>if affine is None:<EOL><INDENT>if isinstance(obj, Subject):<EOL><INDENT>affine = obj.vertex_to_voxel_matrix<EOL><DEDENT>else:<EOL><INDENT>affine = ([[-<NUM_LIT:1>,<NUM_LIT:0>,<NUM_LIT:0>],[<NUM_LIT:0>,<NUM_LIT:0>,-<NUM_LIT:1>],[<NUM_LIT:0>,<NUM_LIT:1>,<NUM_LIT:0>]], np.asarray(shape) / <NUM_LIT:2>)<EOL><DEDENT><DEDENT>affine = to_affine(affine, <NUM_LIT:3>)<EOL>hems = (obj.lh, obj.rh) if isinstance(obj, Subject) else (obj,)<EOL>method = '<STR_LIT>' if method is None else method.lower()<EOL>if method in ['<STR_LIT>', '<STR_LIT>', '<STR_LIT>']:<EOL><INDENT>interp = [_vertex_to_voxel_linear_interpolation(h, mask, shape, affine) for h in hems]<EOL>if len(interp) == <NUM_LIT:1>: interp = interp[<NUM_LIT:0>]<EOL>else: interp = sps.hstack(interp)<EOL><DEDENT>elif method in ['<STR_LIT>', '<STR_LIT>']:<EOL><INDENT>interp = [_vertex_to_voxel_lines_interpolation(h, mask, shape, affine) for h in hems]<EOL>if len(interp) == <NUM_LIT:1>: interp = interp[<NUM_LIT:0>]<EOL>else: interp = sps.hstack(interp)<EOL><DEDENT>elif method in ['<STR_LIT>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>']:<EOL><INDENT>interp = [_vertex_to_voxel_linear_interpolation(h, mask, shape, affine) for h in hems]<EOL>if len(interp) == <NUM_LIT:1>: interp = interp[<NUM_LIT:0>]<EOL>else: interp = sps.hstack(interp)<EOL>(rs,cs) = interp.shape<EOL>argmaxs = np.asarray(interp.argmax(axis=<NUM_LIT:1>))[:,<NUM_LIT:0>]<EOL>return sps.csr_matrix((np.ones(rs, dtype=np.int), (range(rs), argmaxs)),<EOL>shape=interp.shape,<EOL>dtype=np.int)<EOL><DEDENT>elif method in ['<STR_LIT>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>']:<EOL><INDENT>aff = npla.inv(affine)<EOL>interp = [_vertex_to_voxel_nearest_interpolation(h, mask, aff) for h in hems]<EOL>if len(interp) == <NUM_LIT:1>: interp = interp[<NUM_LIT:0>]<EOL>else: interp = sps.hstack(interp)<EOL><DEDENT>else:<EOL><INDENT>raise ValueError('<STR_LIT>' % method)<EOL><DEDENT>return (mask, interp)<EOL>
cortex_to_image_interpolation(obj) yields a tuple (indices, interp) where indices is a tuple of
  voxel indices and interp is an interpolation matrix that converts a vector of cortical
  surface vertex values into a vector of of voxel values with the same ordering as those given
  in indices.

The argument obj may be either a subject or a cortex or a tuple of (lh, rh) cortices. If the obj
specifies two cortices, then the interp matrix will be arranged such that the vertex-value
vector with which interp should be multiplied must first list the LH values then the RH values;
i.e., image[indices] = interp.dot(join(lh_values, rh_values)).

The following options are accepted:
  * mask (default: None) specifies the mask of vertices that should be interpolated; if this is
    None, then this will attempt to use the gray_indices of the object if the object is a
    subject, otherwise this is equivalent to 'all'. The special value 'all' indicates that all
    vertices overlapping with the cortex should be interpolated. Mask may be given as a boolean
    mask or as a set or as a tuple/matrix of indices equivalent to numpy.where(binary_mask) or
    its transpose. The value indices that is returned will be identical to mask if mask is given
    in the same format.
  * affine (default: None) specifies the affine transform that is used to align the vertices
    with the voxels. In voxel-space, the voxel with index (i,j,k) is centered at at the point
    (i,j,k) and is 1 unit wide in every direction. If the value None is given, will attempt to
    use the object's vertex_to_voxel_matrix if object is a subject, otherwise will use a
    FreeSurfer-like orientation that places the vertex origin in the center of the image.
  * shape (default: None) specifies the size of the resulting image as a tuple. If None is given
    then the function attempts to deduce the correct shape; if obj is a subject, then its
    image_dimensions tuple is used. Otherwise, the size is deduced from the mask, if possible;
    if it cannot be deduced from the mask, then (256, 256, 256) is used.
  * method (default: 'linear') specifies the method to use. May be 'linear', 'heaviest', or
    'nearest'.