def __load_images(label_image_n, mask_image_n, boundary_image_n, fg_image_n = False, bg_image_n = False):
    """
    Load and return all image data in preprocessed ndarrays.
    The label image will be relabeled to start from 1.
    @return label, ground-truth-mask, boundary-result-mask[, fg-markers, bg-markers]
    """
    # load images
    label_image = load(label_image_n)
    mask_image = load(mask_image_n)
    boundary_image = load(boundary_image_n)
    if fg_image_n: fg_image = load(fg_image_n)
    if bg_image_n: bg_image = load(bg_image_n)
    
    # extract image data
    label_image_d = scipy.squeeze(label_image.get_data())
    mask_image_d = scipy.squeeze(mask_image.get_data()).astype(scipy.bool_)
    boundary_image_d = scipy.squeeze(boundary_image.get_data()).astype(scipy.bool_)
    if fg_image_n: fg_image_d = scipy.squeeze(fg_image.get_data()).astype(scipy.bool_)
    if bg_image_n: bg_image_d = scipy.squeeze(bg_image.get_data()).astype(scipy.bool_)
    
    # check if images are of same dimensionality
    if label_image_d.shape != mask_image_d.shape: raise argparse.ArgumentError('The mask image {} must be of the same dimensionality as the label image {}.'.format(mask_image_d.shape, label_image_d.shape))
    if label_image_d.shape != boundary_image_d.shape: raise argparse.ArgumentError('The boundary term image {} must be of the same dimensionality as the label image {}.'.format(boundary_image_d.shape, label_image_d.shape))
    if fg_image_n:
        if label_image_d.shape != fg_image_d.shape: raise argparse.ArgumentError('The foreground markers image {} must be of the same dimensionality as the label image {}.'.format(fg_image_d.shape, label_image_d.shape))
    if bg_image_n:
        if label_image_d.shape != bg_image_d.shape: raise argparse.ArgumentError('The background markers image {} must be of the same dimensionality as the label image {}.'.format(bg_image_d.shape, label_image_d.shape))
    
    # relabel the label image to start from 1
    label_image_d = medpy.filter.relabel(label_image_d, 1)
    
    if fg_image_n:
        return label_image_d, mask_image_d, boundary_image_d, fg_image_d, bg_image_d
    else:
        return label_image_d, mask_image_d, boundary_image_d

def test_load_dataarray3():
    img3 = load(DATA_FILE3)
    with InTemporaryDirectory():
        save(img3, "test.gii")
        bimg = load("test.gii")
    for img in (img3, bimg):
        assert_array_almost_equal(img.darrays[0].data[30:50], DATA_FILE3_darr1)

def test_base64_written():
    with InTemporaryDirectory():
        with open(DATA_FILE5, "rb") as fobj:
            contents = fobj.read()
        # Confirm the bad tags are still in the file
        assert_true(b"GIFTI_ENCODING_B64BIN" in contents)
        assert_true(b"GIFTI_ENDIAN_LITTLE" in contents)
        # The good ones are missing
        assert_false(b"Base64Binary" in contents)
        assert_false(b"LittleEndian" in contents)
        # Round trip
        img5 = load(DATA_FILE5)
        save(img5, "fixed.gii")
        with open("fixed.gii", "rb") as fobj:
            contents = fobj.read()
        # The bad codes have gone, replaced by the good ones
        assert_false(b"GIFTI_ENCODING_B64BIN" in contents)
        assert_false(b"GIFTI_ENDIAN_LITTLE" in contents)
        assert_true(b"Base64Binary" in contents)
        if sys.byteorder == "little":
            assert_true(b"LittleEndian" in contents)
        else:
            assert_true(b"BigEndian" in contents)
        img5_fixed = load("fixed.gii")
        darrays = img5_fixed.darrays
        assert_array_almost_equal(darrays[0].data, DATA_FILE5_darr1)
        assert_array_almost_equal(darrays[1].data, DATA_FILE5_darr2)

def test_readwritedata():
    img = load(DATA_FILE2)
    with InTemporaryDirectory():
        save(img, "test.gii")
        img2 = load("test.gii")
        assert_equal(img.numDA, img2.numDA)
        assert_array_almost_equal(img.darrays[0].data, img2.darrays[0].data)

def test_load_dataarray4():
    img4 = load(DATA_FILE4)
    # Round trip
    with InTemporaryDirectory():
        save(img4, "test.gii")
        bimg = load("test.gii")
    for img in (img4, bimg):
        assert_array_almost_equal(img.darrays[0].data[:10], DATA_FILE4_darr1)

def test_load_dataarray2():
    img2 = load(DATA_FILE2)
    # Round trip
    with InTemporaryDirectory():
        save(img2, 'test.gii')
        bimg = load('test.gii')
    for img in (img2, bimg):
        assert_array_almost_equal(img.darrays[0].data[:10], DATA_FILE2_darr1)

def test_read_ordering():
    # DATA_FILE1 has an expected darray[0].data shape of (3,3).  However if we
    # read another image first (DATA_FILE2) then the shape is wrong
    # Read an image
    img2 = load(DATA_FILE2)
    assert_equal(img2.darrays[0].data.shape, (143479, 1))
    # Read image for which we know output shape
    img = load(DATA_FILE1)
    assert_equal(img.darrays[0].data.shape, (3, 3))

def  __load(picklefile, label):
    """
    Load a pickled testbed as well as the original and label image for further processing.
    The label image will be relabeled to start with region id 1.
    @param picklefile the testbed pickle file name
    @param label the label image file name
    @return a tuple containing:
        label: the label image data as ndarray
        bounding_boxes: the bounding boxes around the label image regions (Note that the the bounding box of a region with id rid is accessed using bounding_boxes[rid - 1])
        model_fg_ids: the region ids of all regions to create the foreground model from
        model_bg_ids: the region ids of all regions to create the background model from
        eval_ids: the regions to evaluate the regions term on, represented by their ids
        truth_fg: subset of regions from the eval_ids that are foreground according to the ground-truth
        truth_bg:  subset of regions from the eval_ids that are background according to the ground-truth
    """
    # load and preprocess images
    label_image = load(label)
    
    label_image_d = scipy.squeeze(label_image.get_data())
    
    # relabel the label image to start from 1
    label_image_d = medpy.filter.relabel(label_image_d, 1)
    
    # extracting bounding boxes
    bounding_boxes = find_objects(label_image_d)
    
    # load testbed
    with open(picklefile, 'r') as f:
        model_fg_ids = cPickle.load(f)
        model_bg_ids = cPickle.load(f)
        cPickle.load(f) # eval ids
        truth_fg = cPickle.load(f)
        truth_bg = cPickle.load(f)
            
    return label_image_d, label_image, bounding_boxes, model_fg_ids, model_bg_ids, truth_fg, truth_bg

def test_load_dataarray1():
    img1 = load(DATA_FILE1)
    # Round trip
    with InTemporaryDirectory():
        save(img1, "test.gii")
        bimg = load("test.gii")
    for img in (img1, bimg):
        assert_array_almost_equal(img.darrays[0].data, DATA_FILE1_darr1)
        assert_array_almost_equal(img.darrays[1].data, DATA_FILE1_darr2)
        me = img.darrays[0].meta.metadata
        assert_true("AnatomicalStructurePrimary" in me)
        assert_true("AnatomicalStructureSecondary" in me)
        assert_equal(me["AnatomicalStructurePrimary"], "CortexLeft")
        assert_array_almost_equal(img.darrays[0].coordsys.xform, np.eye(4, 4))
        assert_equal(xform_codes.niistring[img.darrays[0].coordsys.dataspace], "NIFTI_XFORM_TALAIRACH")
        assert_equal(xform_codes.niistring[img.darrays[0].coordsys.xformspace], "NIFTI_XFORM_TALAIRACH")

def test_load_labeltable():
    img6 = load(DATA_FILE6)
    # Round trip
    with InTemporaryDirectory():
        save(img6, "test.gii")
        bimg = load("test.gii")
    for img in (img6, bimg):
        assert_array_almost_equal(img.darrays[0].data[:3], DATA_FILE6_darr1)
        assert_equal(len(img.labeltable.labels), 36)
        labeldict = img.labeltable.get_labels_as_dict()
        assert_true(660700 in labeldict)
        assert_equal(labeldict[660700], "entorhinal")
        assert_equal(img.labeltable.labels[1].key, 2647065)
        assert_equal(img.labeltable.labels[1].red, 0.0980392)
        assert_equal(img.labeltable.labels[1].green, 0.392157)
        assert_equal(img.labeltable.labels[1].blue, 0.156863)
        assert_equal(img.labeltable.labels[1].alpha, 1)

def test_save_load():
    shape = (2, 4, 6)
    npt = np.float32
    data = np.arange(np.prod(shape), dtype=npt).reshape(shape)
    affine = np.diag([1, 2, 3, 1])
    affine[:3,3] = [3,2,1]
    img = ni1.Nifti1Image(data, affine)
    img.set_data_dtype(npt)
    with InTemporaryDirectory() as pth:
        pth = mkdtemp()
        nifn = pjoin(pth, 'an_image.nii')
        sifn = pjoin(pth, 'another_image.img')
        ni1.save(img, nifn)
        re_img = nils.load(nifn)
        yield assert_true(isinstance(re_img, ni1.Nifti1Image))
        yield assert_array_equal(re_img.get_data(), data)
        yield assert_array_equal(re_img.get_affine(), affine)
        # These and subsequent del statements are to prevent confusing
        # windows errors when trying to open files or delete the
        # temporary directory. 
        del re_img
        try:
            import scipy.io
        except ImportError:
            # ignore if there is no matfile reader, and restart
            pass
        else:
            spm2.save(img, sifn)
            re_img2 = nils.load(sifn)
            yield assert_true(isinstance(re_img2, spm2.Spm2AnalyzeImage))
            yield assert_array_equal(re_img2.get_data(), data)
            yield assert_array_equal(re_img2.get_affine(), affine)
            del re_img2
            spm99.save(img, sifn)
            re_img3 = nils.load(sifn)
            yield assert_true(isinstance(re_img3,
                                         spm99.Spm99AnalyzeImage))
            yield assert_array_equal(re_img3.get_data(), data)
            yield assert_array_equal(re_img3.get_affine(), affine)
            ni1.save(re_img3, nifn)
            del re_img3
        re_img = nils.load(nifn)
        yield assert_true(isinstance(re_img, ni1.Nifti1Image))
        yield assert_array_equal(re_img.get_data(), data)
        yield assert_array_equal(re_img.get_affine(), affine)
        del re_img

def test_parse_dataarrays():
    fn = "bad_daa.gii"
    img = gi.GiftiImage()

    with InTemporaryDirectory():
        save(img, fn)
        with open(fn, "r") as fp:
            txt = fp.read()
        # Make a bad gifti.
        txt = txt.replace('NumberOfDataArrays="0"', 'NumberOfDataArrays ="1"')
        with open(fn, "w") as fp:
            fp.write(txt)

        with clear_and_catch_warnings() as w:
            warnings.filterwarnings("once", category=UserWarning)
            load(fn)
            assert_equal(len(w), 1)
            assert_equal(img.numDA, 0)

def main():
    # parse cmd arguments
    parser = getParser()
    parser.parse_args()
    args = getArguments(parser)
    
    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)
    
    # build output image name
    image_bg_name = args.folder + '/' + args.mask.split('/')[-1][:-4] + '.bg'
    image_bg_name += args.mask.split('/')[-1][-4:]
        
    # check if output image exists
    if not args.force:
        if os.path.exists(image_bg_name):
            logger.warning('The output image {} already exists. Breaking.'.format(image_bg_name))
            exit(1)
    
    # load mask
    logger.info('Loading mask {}...'.format(args.mask))
    try: 
        mask_image = load(args.mask)
        mask_image_data = numpy.squeeze(mask_image.get_data()).astype(scipy.bool_)
    except ImageFileError as e:
        logger.critical('The mask image does not exist or its file type is unknown.')
        raise ArgumentError('The mask image does not exist or its file type is unknown.', e)  
    
    # array of indices to access desired slices
    sls = [(slice(1), slice(None), slice(None)),
           (slice(-1, None), slice(None), slice(None)),
           (slice(None), slice(1), slice(None)),
           (slice(None), slice(-1, None), slice(None)),
           (slice(None), slice(None), slice(1)),
           (slice(None), slice(None), slice(-1, None))]
    
    # security check
    logger.info('Determine if the slices are not intersection with the reference liver mask...')
    for sl in sls:
        if not 0 == len(mask_image_data[sl].nonzero()[0]):
            logger.critical('Reference mask reaches till the image border.')
            raise ArgumentError('Reference mask reaches till the image border.')
        
    # create and save background marker image
    logger.info('Creating background marker image...')
    image_bg_data = scipy.zeros(mask_image_data.shape, dtype=scipy.bool_)
    for sl in sls:
        image_bg_data[sl] = True
    
    logger.info('Saving background marker image...')
    mask_image.get_header().set_data_dtype(scipy.int8)
    save(image_like(image_bg_data, mask_image), image_bg_name)
    
    logger.info('Successfully terminated.')

def test_save_load():
    shape = (2, 4, 6)
    npt = np.float32
    data = np.arange(np.prod(shape), dtype=npt).reshape(shape)
    affine = np.diag([1, 2, 3, 1])
    affine[:3,3] = [3,2,1]
    img = ni1.Nifti1Image(data, affine)
    img.set_data_dtype(npt)
    with InTemporaryDirectory() as pth:
        nifn = 'an_image.nii'
        sifn = 'another_image.img'
        ni1.save(img, nifn)
        re_img = nils.load(nifn)
        yield assert_true(isinstance(re_img, ni1.Nifti1Image))
        yield assert_array_equal(re_img.get_data(), data)
        yield assert_array_equal(re_img.get_affine(), affine)
        # These and subsequent del statements are to prevent confusing
        # windows errors when trying to open files or delete the
        # temporary directory. 
        del re_img
        if have_scipy: # skip we we cannot read .mat files
            spm2.save(img, sifn)
            re_img2 = nils.load(sifn)
            yield assert_true(isinstance(re_img2, spm2.Spm2AnalyzeImage))
            yield assert_array_equal(re_img2.get_data(), data)
            yield assert_array_equal(re_img2.get_affine(), affine)
            del re_img2
            spm99.save(img, sifn)
            re_img3 = nils.load(sifn)
            yield assert_true(isinstance(re_img3,
                                         spm99.Spm99AnalyzeImage))
            yield assert_array_equal(re_img3.get_data(), data)
            yield assert_array_equal(re_img3.get_affine(), affine)
            ni1.save(re_img3, nifn)
            del re_img3
        re_img = nils.load(nifn)
        yield assert_true(isinstance(re_img, ni1.Nifti1Image))
        yield assert_array_equal(re_img.get_data(), data)
        yield assert_array_equal(re_img.get_affine(), affine)
        del re_img

def test_metadata_deprecations():
    img = load(datafiles[0])
    me = img.meta

    # Test deprecation
    with clear_and_catch_warnings() as w:
        warnings.filterwarnings("once", category=DeprecationWarning)
        assert_equal(me, img.get_meta())

    with clear_and_catch_warnings() as w:
        warnings.filterwarnings("once", category=DeprecationWarning)
        img.set_metadata(me)
    assert_equal(me, img.meta)

def test_labeltable_deprecations():
    img = load(DATA_FILE6)
    lt = img.labeltable

    # Test deprecation
    with clear_and_catch_warnings() as w:
        warnings.filterwarnings('once', category=DeprecationWarning)
        assert_equal(lt, img.get_labeltable())

    with clear_and_catch_warnings() as w:
        warnings.filterwarnings('once', category=DeprecationWarning)
        img.set_labeltable(lt)
    assert_equal(lt, img.labeltable)

def main():
    # prepare logger
    a = [[[1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5]],
         [[1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5]],
         [[1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5],
          [1,2,3,4,5]]]
    
    b = [[[1,2,3,4,3],
          [3,1,2,1,4],
          [3,1,4,2,2],
          [2,2,1,3,3],
          [4,3,1,2,1]]]
    
    c = [[[1,2,3,4],
          [5,6,7,8],
          [9,1,2,3]],
         [[4,5,6,7],
          [8,9,1,2],
          [3,4,5,6]]]
    
    a = scipy.asarray(a)
    b = scipy.asarray(b)
    c = scipy.asarray(c)
    
    original_image_n = '/home/omaier/Experiments/Regionsegmentation/Evaluation_Viscous/00originalvolumes/o09.nii'
    original_image = load(original_image_n)
    original_image_d = scipy.squeeze(original_image.get_data())
    
    original_image_d = original_image_d[0:50, 0:50, 0:50]
    
    dir = scipy.zeros(original_image_d.shape).ravel()
    oir = original_image_d.ravel()
    for i in range(len(dir)):
        dir[i] = oir[i]
    dir = dir.reshape(original_image_d.shape)
    
    coa, con, dir = tamura(dir, 5)
    
    for i in range(5):
        print "k=", i+1, " with ", len((coa == i).nonzero()[0])

def main():
    # parse cmd arguments
    parser = getParser()
    parser.parse_args()
    args = getArguments(parser)
    
    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)
    
    # build output file name
    file_csv_name = args.csv + '.csv'
    
    # check if output file exists
    if not args.force:
        if os.path.exists(file_csv_name):
            logger.warning('The output file {} already exists. Skipping.'.format(file_csv_name))
            sys.exit(0)
    
    # open output file
    with open(file_csv_name, 'w') as f:
        
        # write header into file
        f.write('image;min_x;min_y;min_z;max_x;max_y;max_z\n')
        
        # iterate over input images
        for image in args.images:
            
            # get and prepare image data
            logger.info('Processing image {}...'.format(image))
            image_data = numpy.squeeze(load(image).get_data())
            
            # count number of labels and flag a warning if they reach the ushort border
            mask = image_data.nonzero()

            # count number of labels and write into file
            f.write('{};{};{};{};{};{};{}\n'.format(image.split('/')[-1],
                                                    mask[0].min(),
                                                    mask[1].min(),
                                                    mask[2].min(),
                                                    mask[0].max(),
                                                    mask[1].max(),
                                                    mask[2].max()))
            
            f.flush()
            
    logger.info('Successfully terminated.')

def test_load_getbyintent():
    img = load(DATA_FILE1)
    da = img.get_arrays_from_intent("NIFTI_INTENT_POINTSET")
    assert_equal(len(da), 1)

    with clear_and_catch_warnings() as w:
        warnings.filterwarnings("once", category=DeprecationWarning)
        da = img.getArraysFromIntent("NIFTI_INTENT_POINTSET")
        assert_equal(len(da), 1)
        assert_equal(len(w), 1)
        assert_equal(w[0].category, DeprecationWarning)

    da = img.get_arrays_from_intent("NIFTI_INTENT_TRIANGLE")
    assert_equal(len(da), 1)

    da = img.get_arrays_from_intent("NIFTI_INTENT_CORREL")
    assert_equal(len(da), 0)
    assert_equal(da, [])

def main():
    # prepare logger
    logger = Logger.getInstance()
    logger.setLevel(logging.DEBUG)
    
    # input image locations
    #i = '/home/omaier/Experiments/Regionsegmentation/Evaluation_Viscous/00originalvolumes/o09.nii' # original image
    #i = '/home/omaier/Temp/test.nii' # original image
    #i = '/home/omaier/Temp/o09_smoothed_i4.0_c0.1_t0.0625.nii'
    i = '/home/omaier/Experiments/GraphCut/BoundaryTerm/Stawiaski/01gradient/o09_gradient.nii'
    
    # output image locations
    r = '/home/omaier/Temp/result_gradient.nii' # result mask
    
    # load images
    i_i = load(i)
 
    # extract and prepare image data
    i_d = scipy.squeeze(i_i.get_data())
    
    # crop input images to achieve faster execution
    crop = [slice(50, -200),
            slice(50, -150),
            slice(50, -100)]
    #i_d = i_d[crop]   
    
    i_d = scipy.copy(i_d)
    
    # !TODO: Test if input image is of size 0
    
    logger.debug('input image shape={},ndims={},dtype={}'.format(i_d.shape, i_d.ndim, i_d.dtype))

    result = watershed8(i_d, logger)

    logger.info('Saving resulting region map...')
    result_i = image_like(result, i_i)
    result_i.get_header().set_data_dtype(scipy.int32)
    save(result_i, r)

    logger.info('Done!')