docs: def _read_image_data( mhd_file):
      """ Read the resulting image data and return it as a numpy array.
      """
      tempdir = get_tempdir()
      
      # Load description from mhd file
      fname = tempdir + '/' + mhd_file
      des = open(fname, 'r').read()
      
      # Get data filename and load raw data
      match = re.findall('ElementDataFile = (.+?)\n', des)
      fname = tempdir + '/' + match[0]
      data = open(fname, 'rb').read()
      
      # Determine dtype
      match = re.findall('ElementType = (.+?)\n', des)
      dtype_itk = match[0].upper().strip()
      dtype = DTYPE_ITK2NP.get(dtype_itk, None)
      if dtype is None:
          raise RuntimeError('Unknown ElementType: ' + dtype_itk)
      
      # Create numpy array
      a = np.frombuffer(data, dtype=dtype)
      
      # Determine shape, sampling and origin of the data
      match = re.findall('DimSize = (.+?)\n', des)
      shape = [int(i) for i in match[0].split(' ')]
      #
      match = re.findall('ElementSpacing = (.+?)\n', des)
      sampling = [float(i) for i in match[0].split(' ')]
      #
      match = re.findall('Offset = (.+?)\n', des)
      origin = [float(i) for i in match[0].split(' ')]
      
      # Reverse shape stuff to make z-y-x order
      shape = [s for s in reversed(shape)]
      sampling = [s for s in reversed(sampling)]
      origin = [s for s in reversed(origin)]
      
      # Take vectors/colours into account
      N = np.prod(shape)
      if N != a.size:
          extraDim = int( a.size / N )
          shape = tuple(shape) + (extraDim,)
          sampling = tuple(sampling) + (1.0,)
          origin = tuple(origin) + (0,)
      
      # Check shape
      N = np.prod(shape)
      if N != a.size:
          raise RuntimeError('Cannot apply shape to data.')
      else:
          a.shape = shape
          a = Image(a)
          a.sampling = sampling
          a.origin = origin
      return a