def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)

    # load input images
    input_data, input_header = load(args.input)
    original_data, _ = load(args.original)
    
    logger.debug('Old shape={}.'.format(input_data.shape))
    
    # compute position
    logger.info('Computing positon and pad volume...')
    position = __parse_contour_list(args.contours, input_data)
    
    # pad volume
    output_data = scipy.zeros(original_data.shape, input_data.dtype)
    output_data[position] = input_data
    
    
    logger.debug('New shape={}.'.format(input_data.shape))
    
    # save result contour volume
    save(output_data, args.output, input_header, args.force)

    logger.info("Successfully terminated.")

def sresamplebyexample(src, dest, referenceimage, binary = False):
    r"""
    Secure-re-sample an image located at ``src`` by example ``referenceimage`` and
    save it under ``dest``.
    
    Parameters
    ----------
    src : string
        Source image file.
    dest : string
        Destination image file.
    referenceimage : string
        Reference image displaying the target spacing, origin and size.
    binary : bool
        Set to ``True`` for binary images.
    """
    # get target voxel spacing
    refimage, refhdr = load(referenceimage)
    spacing = header.get_pixel_spacing(refhdr)
    
    with tmpdir() as t:
        # create a temporary copy of the reference image with the source image data-type (imiImageResample requires both images to be of the same dtype)
        srcimage, _ = load(src)
        save(refimage.astype(srcimage.dtype), os.path.join(t, 'ref.nii.gz'), refhdr)
    
        # prepare and run registration command
        cmd = ['imiImageResample', '-I', src, '-O', dest, '-R', os.path.join(t, 'ref.nii.gz'), '-s'] + map(str, spacing)
        if binary:
            cmd += ['-b']
        rtcode, stdout, stderr = call(cmd)
    
    # check if successful
    if not os.path.isfile(dest):
        raise CommandExecutionError(cmd, rtcode, stdout, stderr, 'Binary re-sampling result image not created.')

def main():
	i1, h1 = load(sys.argv[1])
	i2, h2 = load(sys.argv[2])

	# shift image to align origins
	origin_h1 = numpy.sign(h1.get_qform()[0:3,0:3]).dot(header.get_offset(h1))
	origin_h2 = numpy.sign(h2.get_qform()[0:3,0:3]).dot(header.get_offset(h2))
	origin_difference_pixel = (origin_h1 - origin_h2) / numpy.asarray(header.get_pixel_spacing(h1))
	# negative values: shift image 1 by this upon inserting (which is the smae as cutting the output image)
	# positive values: cut image 1 by this at inserting and also cut right side by length of output image plus this value
	o = numpy.zeros(i2.shape, i2.dtype)
	o_slicer = []
	i_slicer = []
	for j, p in enumerate(origin_difference_pixel):
		if p >= 0:
			i_slicer.append(slice(0,      min(i1.shape[j], o.shape[j] - abs(p))))
			o_slicer.append(slice(abs(p), min(i1.shape[j] + abs(p), o.shape[j])))
		else:
			i_slicer.append(slice(abs(p), min(i1.shape[j], o.shape[j] + abs(p))))
			o_slicer.append(slice(0,      min(i1.shape[j] - abs(p), o.shape[j])))

	o[o_slicer] = i1[i_slicer]
	header.set_offset(h1, header.get_offset(h2))
	
	save(o, sys.argv[3], h1)

def main():
	onedir = sys.argv[1] # the first folder containing case folders
	twodir = sys.argv[2] # the second folder containing case folders
	nocase = (len(sys.argv) > 3 and sys.argv[3] == '-i')

	if DEBUG: print 'INFO: Comparing all cases in folders {} and {}.'.format(onedir, twodir)

	# iterate over first folder and compare voxel spacings with equivalent image in second folder
	print "Case\tvs same\tshape same"
	for root, dirs, files in os.walk(onedir):
		for case in sorted(dirs):
			for root, dirs, files in os.walk('{}/{}'.format(onedir, case)):
				for file_ in files:
					if file_.endswith(FILE_ENDING):
						i, hi = load('{}/{}/{}'.format(onedir, case, file_))
						if nocase:
							j, hj = load('{}/{}.{}'.format(twodir, case, FILE_ENDING))
						else:
							j, hj = load('{}/{}/{}'.format(twodir, case, file_))
						vs_same = numpy.array_equal(header.get_pixel_spacing(hi), header.get_pixel_spacing(hj))
						shape_same = numpy.array_equal(i.shape, j.shape)
						print '{}\t{}\t{}'.format(case, vs_same, shape_same)
						if not vs_same:
							print "\t{} vs {}".format(header.get_pixel_spacing(hi), header.get_pixel_spacing(hj))
						if not shape_same:
							print "\t{} vs {}".format(i.shape, j.shape)
	print 'Terminated.'

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)
    
    # load input image1
    data_input1, _ = load(args.input1)
    
    # load input image2
    data_input2, _ = load(args.input2)
    
    # compare dtype and shape
    if not data_input1.dtype == data_input2.dtype: print 'Dtype differs: {} to {}'.format(data_input1.dtype, data_input2.dtype)
    if not data_input1.shape == data_input2.shape:
        print 'Shape differs: {} to {}'.format(data_input1.shape, data_input2.shape)
        print 'The voxel content of images of different shape can not be compared. Exiting.'
        sys.exit(-1)
    
    # compare image data
    voxel_total = reduce(lambda x, y: x*y, data_input1.shape)
    voxel_difference = len((data_input1 != data_input2).nonzero()[0])
    if not 0 == voxel_difference:
        print 'Voxel differ: {} of {} total voxels'.format(voxel_difference, voxel_total)
        print 'Max difference: {}'.format(scipy.absolute(data_input1 - data_input2).max())
    else: print 'No other difference.'
    
    logger.info("Successfully terminated.")    

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)  
    
    # constants
    colours = {'i': 10, 'o': 11}
    
    # load volumes
    marker_data, _ = load(args.marker)
    contour_data, _ = load(args.contour)
    
    # perform check
    contour_data = contour_data == colours[args.type]
    marker_data_fg = marker_data == 1
    marker_data_bg = marker_data == 2
    if scipy.logical_and(contour_data, marker_data_fg).any():
        logger.warning('Intersection between {} and {} (type {}) in foreground.'.format(args.marker, args.contour, args.type))
    elif scipy.logical_and(contour_data, marker_data_bg).any():
        logger.warning('Intersection between {} and {} (type {}) in background.'.format(args.marker, args.contour, args.type))
    else:
        print "No intersection."

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)
    
    # load input image using nibabel
    logger.info('Loading image {}...'.format(args.input))
    image_labels_data, _ = load(args.image)    
    
    # load mask image
    logger.info('Loading mask {}...'.format(args.mask))
    image_mask_data, image_mask_data_header = load(args.mask)
    
    # check if output image exists
    if not args.force:
        if os.path.exists(args.output):
            logger.warning('The output image {} already exists. Skipping this image.'.format(args.output))
    
    # create a mask from the label image
    logger.info('Reducing the label image...')
    image_reduced_data = fit_labels_to_mask(image_labels_data, image_mask_data)
    
    # save resulting mask
    logger.info('Saving resulting mask as {} in the same format as input mask, only with data-type int8...'.format(args.output))
    image_reduced_data = image_reduced_data.astype(numpy.bool, copy=False) # bool sadly not recognized
    save(image_reduced_data, args.output, image_mask_data_header, args.force)
    
    logger.info('Successfully terminated.')

def test03(img, hdr, idx, delta):
    # TEST 03: DOES ANY META-INFORMATION GET LOST DURING FORMAT CONVERSION? AND IF YES; WHICH?
    for tr in types_int: # reference type
        print ''
        oformat = tr
        
        # create, save and load reference image
        try:
            name_ref = tmp_folder + '.'.join(['tmp_ref', tr])
            save(img, name_ref, hdr, True)
            img_ref, hdr_ref = load(name_ref)
        except Exception as e:
            print '\tERROR: Could not generate reference image for type {}: {}'.format(otype, e)
            continue
        
        # extract meta-data from reference image
        mdata_ref = {'shape': img_ref.shape,
                     'dtype': img_ref.dtype,
                     'point': img_ref[idx],
                     'spacing': header.get_pixel_spacing(hdr_ref),
                     'offset': header.get_offset(hdr_ref),}        
        
        # print meta-data from reference image
        
        # iterate of test images
        for tt in types_int: # test type
            print '{} => {}'.format(oformat, tt),
            
            # create, save and load test images
            try:
                #print type(img_ref), type(hdr_ref)
                #print type(img_test), type(hdr_test)
                name_test = tmp_folder + '.'.join(['tmp_test', tt])
                save(img_ref, name_test, hdr_ref, True)
                img_test, hdr_test = load(name_test)
                
            except Exception as e:
                print '\tERROR: Could not generate test image. {}'.format(e)
                continue
            
            # extract meta-data from test image
            mdata_test = {'shape': img_test.shape,
                          'dtype': img_test.dtype,
                          'spacing': header.get_pixel_spacing(hdr_test),
                          'offset': header.get_offset(hdr_test),
                          'point': img_test[idx]}                    
            
            # compare reference against meta-image
            error = False
            for k in mdata_ref.keys():
                equal = _compare(mdata_ref[k], mdata_test[k], delta)
                #print '\n\t{} ({}) : {} = {}'.format(equal, k, mdata_ref[k], mdata_test[k]),
                if not equal:
                    error = True
                    print '\n\t{} ({}) : {} = {}'.format(equal, k, mdata_ref[k], mdata_test[k]),
            if not error:
                print '\t{}'.format(True)
            else:
                print '\n'

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)
    
    # loading input images
    b0img, b0hdr = load(args.b0image)
    bximg, bxhdr = load(args.bximage)
    
    # convert to float
    b0img = b0img.astype(numpy.float)
    bximg = bximg.astype(numpy.float)

    # check if image are compatible
    if not b0img.shape == bximg.shape:
        raise ArgumentError('The input images shapes differ i.e. {} != {}.'.format(b0img.shape, bximg.shape))
    if not header.get_pixel_spacing(b0hdr) == header.get_pixel_spacing(bxhdr):
        raise ArgumentError('The input images voxel spacing differs i.e. {} != {}.'.format(header.get_pixel_spacing(b0hdr), header.get_pixel_spacing(bxhdr)))
    
    # check if supplied threshold value as well as the b value is above 0
    if args.threshold is not None and not args.threshold >= 0:
        raise ArgumentError('The supplied threshold value must be greater than 0, otherwise a division through 0 might occur.')
    if not args.b > 0:
        raise ArgumentError('The supplied b-value must be greater than 0.')
    
    # compute threshold value if not supplied
    if args.threshold is None:
        b0thr = otsu(b0img, 32) / 4. # divide by 4 to decrease impact
        bxthr = otsu(bximg, 32) / 4.
        if 0 >= b0thr:
            raise ArgumentError('The supplied b0image seems to contain negative values.')
        if 0 >= bxthr:
            raise ArgumentError('The supplied bximage seems to contain negative values.')
    else:
        b0thr = bxthr = args.threshold
    
    logger.debug('thresholds={}/{}, b-value={}'.format(b0thr, bxthr, args.b))
    
    # threshold b0 + bx DW image to obtain a mask
    # b0 mask avoid division through 0, bx mask avoids a zero in the ln(x) computation
    mask = binary_fill_holes(b0img > b0thr) & binary_fill_holes(bximg > bxthr)
    
    # perform a number of binary morphology steps to select the brain only
    mask = binary_erosion(mask, iterations=1)
    mask = largest_connected_component(mask)
    mask = binary_dilation(mask, iterations=1)
    
    logger.debug('excluding {} of {} voxels from the computation and setting them to zero'.format(numpy.count_nonzero(mask), numpy.prod(mask.shape)))
    
    # compute the ADC
    adc = numpy.zeros(b0img.shape, b0img.dtype)
    adc[mask] = -1. * args.b * numpy.log(bximg[mask] / b0img[mask])
    adc[adc < 0] = 0
            
    # saving the resulting image
    save(adc, args.output, b0hdr, args.force)

def main():
	m = load(sys.argv[1])[0].astype(numpy.bool)
	s = load(sys.argv[2])[0].astype(numpy.bool)

	intc = numpy.count_nonzero(~m & s)

	print "Non-intersecting part of the segmentation:"
	print "{} out of {} voxels".format(intc, numpy.count_nonzero(s))

def main():
	# load input image
	i, _ = load(sys.argv[1])

	# load template image
	_, h = load(sys.argv[2])
	
	# save input image with adapted header in place
	j = i.copy()
	save(j, sys.argv[1], h)

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)
    
    # load first input image as example 
    example_data, example_header = load(args.inputs[0])
    
    # test if the supplied position is valid
    if args.position > example_data.ndim or args.position < 0:
        raise ArgumentError('The supplied position for the new dimension is invalid. It has to be between 0 and {}.'.format(example_data.ndim))
    
    # prepare empty output volume
    output_data = scipy.zeros([len(args.inputs)] + list(example_data.shape), dtype=example_data.dtype)
    
    # add first image to output volume
    output_data[0] = example_data
    
    # load input images and add to output volume
    for idx, image in enumerate(args.inputs[1:]):
        image_data, _ = load(image)
        if not args.ignore and image_data.dtype != example_data.dtype:
            raise ArgumentError('The dtype {} of image {} differs from the one of the first image {}, which is {}.'.format(image_data.dtype, image, args.inputs[0], example_data.dtype))
        if image_data.shape != example_data.shape:
            raise ArgumentError('The shape {} of image {} differs from the one of the first image {}, which is {}.'.format(image_data.shape, image, args.inputs[0], example_data.shape))
        output_data[idx + 1] = image_data
        
    # move new dimension to the end or to target position
    for dim in range(output_data.ndim - 1):
        if dim >= args.position: break
        output_data = scipy.swapaxes(output_data, dim, dim + 1)
        
    # set pixel spacing
    spacing = list(header.get_pixel_spacing(example_header))
    spacing = tuple(spacing[:args.position] + [args.spacing] + spacing[args.position:])
    
    # !TODO: Find a way to enable this also for PyDicom and ITK images
    if __is_header_nibabel(example_header):
        __update_header_from_array_nibabel(example_header, output_data)
        header.set_pixel_spacing(example_header, spacing)
    else:
        raise ArgumentError("Sorry. Setting the voxel spacing of the new dimension only works with NIfTI images. See the description of this program for more details.")
    
    # save created volume
    save(output_data, args.output, example_header, args.force)
        
    logger.info("Successfully terminated.")

def main():
	# loading the image mask
	m = load(sys.argv[2])[0].astype(numpy.bool)

	# extracting the required features and saving them
	for sequence, function_call, function_arguments, voxelspacing in features_to_extract:
		if not isfv(sys.argv[3], sequence, function_call, function_arguments):
			#print sequence, function_call.__name__, function_arguments
			i, h = load('{}/{}.nii.gz'.format(sys.argv[1], sequence))
			call_arguments = list(function_arguments)
			if voxelspacing: call_arguments.append(header.get_pixel_spacing(h))
			call_arguments.append(m)
			fv = function_call(i, *call_arguments)
			savefv(fv, sys.argv[3], sequence, function_call, function_arguments)

def _percentilemodelstandardisation(trainingfiles, brainmaskfiles, destfiles, destmodel):
    r"""
    Train an intensity standardisation model and apply it. All values outside of the
    brain mask are set to zero.
    
    Parameters
    ----------
    trainingfiles : sequence of strings
        All images to use for training and to which subsequently apply the trained model.
    brainmaskfiles : sequence of strings
        The brain masks corresponding to ``trainingfiles``.
    destfiles : sequence of strings
        The intensity standarised target locations corresponding to ``trainingfiles``.
    destmodel : string
        The target model location.
    """
    # check arguments
    if not len(trainingfiles) == len(brainmaskfiles):
        raise ValueError('The number of supplied trainingfiles must be equal to the number of brainmaskfiles.')
    elif not len(trainingfiles) == len(destfiles):
        raise ValueError('The number of supplied trainingfiles must be equal to the number of destfiles.')
    
    # loading input images (as image, header pairs)
    images = []
    headers = []
    for image_name in trainingfiles:
        i, h = load(image_name)
        images.append(i)
        headers.append(h)
        
    # loading brainmasks
    masks = [load(mask_name)[0].astype(numpy.bool) for mask_name in brainmaskfiles]
        
    # train the model
    irs = IntensityRangeStandardization()
    trained_model, transformed_images = irs.train_transform([i[m] for i, m in zip(images, masks)])
    
    # condense outliers in the image (extreme peak values at both end-points of the histogram)
    transformed_images = [_condense(i) for i in transformed_images]
    
    # saving the model
    with open(destmodel, 'wb') as f:
        pickle.dump(trained_model, f)
    
    # save the transformed images
    for ti, i, m, h, dest in zip(transformed_images, images, masks, headers, destfiles):
        i[~m] = 0
        i[m] = ti
        save(i, dest, h)

def main():

	# catch parameters
	segmentation_base_string = sys.argv[1]
	ground_truth_base_string = sys.argv[2]
	mask_file_base_string = sys.argv[3]
	cases = sys.argv[4:]

	# evaluate each case and collect the scores
	hds = []
	assds = []
	precisions = []
	recalls = []
	dcs = []

	# load images and apply mask to segmentation and ground truth (to remove ground truth fg outside of brain mask)
	splush = [load(segmentation_base_string.format(case)) for case in cases]
	tplush = [load(ground_truth_base_string.format(case)) for case in cases]
	masks = [load(mask_file_base_string.format(case))[0].astype(numpy.bool) for case in cases]

	s = [s.astype(numpy.bool) & m for (s, _), m in zip(splush, masks)]
	t = [t.astype(numpy.bool) & m for (t, _), m in zip(tplush, masks)]
	hs = [h for _, h in splush]
	ht = [h for _, h in tplush]

	# compute and append metrics (Pool-processing)
	pool = Pool(n_jobs)
	dcs = pool.map(wdc, zip(t, s))
	precisions = pool.map(wprecision, zip(s, t))
	recalls = pool.map(wrecall, zip(s, t))
	hds = pool.map(whd, zip(t, s, [header.get_pixel_spacing(h) for h in ht]))
	assds = pool.map(wassd, zip(t, s, [header.get_pixel_spacing(h) for h in ht]))

	# print case-wise results
	print 'Metrics:'
	print 'Case\tDC[0,1]\tHD(mm)\tP2C(mm)\tprec.\trecall'
    	for case, _dc, _hd, _assd, _pr, _rc in zip(cases, dcs, hds, assds, precisions, recalls):
        	print '{}\t{:>3,.3f}\t{:>4,.3f}\t{:>4,.3f}\t{:>3,.3f}\t{:>3,.3f}'.format(case, _dc, _hd, _assd, _pr, _rc)
        
	# check for nan/inf values of failed cases and signal warning
	mask = numpy.isfinite(hds)
	if not numpy.all(mask):
		print 'WARNING: Average values only computed on {} of {} cases!'.format(numpy.count_nonzero(mask), mask.size)
		
    	print 'DM  average\t{} +/- {} (Median: {})'.format(numpy.asarray(dcs)[mask].mean(), numpy.asarray(dcs)[mask].std(), numpy.median(numpy.asarray(dcs)[mask]))
    	print 'HD  average\t{} +/- {} (Median: {})'.format(numpy.asarray(hds)[mask].mean(), numpy.asarray(hds)[mask].std(), numpy.median(numpy.asarray(hds)[mask]))
    	print 'ASSD average\t{} +/- {} (Median: {})'.format(numpy.asarray(assds)[mask].mean(), numpy.asarray(assds)[mask].std(), numpy.median(numpy.asarray(assds)[mask]))
    	print 'Prec. average\t{} +/- {} (Median: {})'.format(numpy.asarray(precisions)[mask].mean(), numpy.asarray(precisions)[mask].std(), numpy.median(numpy.asarray(precisions)[mask]))
    	print 'Rec. average\t{} +/- {} (Median: {})'.format(numpy.asarray(recalls)[mask].mean(), numpy.asarray(recalls)[mask].std(), numpy.median(numpy.asarray(recalls)[mask]))

def main():
    i, h = load(sys.argv[1])
    r, _ = load(sys.argv[2])

    diff = numpy.asarray(r.shape) - numpy.asarray(i.shape)

    if numpy.any(diff < 0): # cut to fit
        slicers = [slice(None) if 0 == e else slice(abs(e) / 2, -1 * (abs(e) / 2 + abs(e) % 2)) for e in diff]
        o = i[slicers]
        
    else: # pad to fit
        padding = [(e / 2, e / 2 + e % 2) for e in diff]
        o = numpy.pad(i, padding, "constant")

    save(o, sys.argv[1], h, True)

    def _run_interface(self, runtime):
        if not base.isdefined(self.inputs.out_file):
            self.inputs.out_file = self._gen_filename('out_file')

        in_file = self.inputs.in_file
        mask_file = self.inputs.mask_file
        out_file = self.inputs.out_file

        image, header = mio.load(in_file)
        mask, _ = mio.load(mask_file)

        image[~(mask.astype(numpy.bool))] = 0
        mio.save(image, out_file, header)

        return runtime

def main():
	print 'lesion\tvolume (%)\tvolume (mm)'

	files = [f for f in os.listdir('{}'.format(sys.argv[1])) if os.path.isfile('{}/{}'.format(sys.argv[1], f))]
	for f in files:
		l, h = load('{}/{}'.format(sys.argv[1], f))
		m, _ = load('{}/{}'.format(sys.argv[2], f))

		lesion_voxel = numpy.count_nonzero(l)
		total_voxel = numpy.count_nonzero(m)

		volume_mm = numpy.prod(header.get_pixel_spacing(h)) * lesion_voxel
		volume_percentage = lesion_voxel / float(total_voxel)

		print '{}\t{}\t{}\t'.format(f[:-7], volume_percentage, volume_mm)

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug: logger.setLevel(logging.DEBUG)
    elif args.verbose: logger.setLevel(logging.INFO)
    
    # loading input images (as image, header pairs)
    images = []
    headers = []
    for image_name in args.images:
        i, h = load(image_name)
        images.append(i)
        headers.append(h)
    
    # loading binary foreground masks if supplied, else create masks from threshold value
    if args.masks:
        masks = [load(mask_name)[0].astype(numpy.bool) for mask_name in args.masks]
    else:
        masks = [i > args.threshold for i in images]
    
    # if in application mode, load the supplied model and apply it to the images
    if args.lmodel:
        logger.info('Loading the model and transforming images...')
        with open(args.lmodel, 'r') as f:
            trained_model = pickle.load(f)
            if not isinstance(trained_model, IntensityRangeStandardization):
                raise ArgumentError('{} does not seem to be a valid pickled instance of an IntensityRangeStandardization object'.format(args.lmodel))
            transformed_images = [trained_model.transform(i[m], surpress_mapping_check = args.ignore) for i, m in zip(images, masks)]
            
    # in in training mode, train the model, apply it to the images and save it
    else:
        logger.info('Training the average intensity model...')
        irs = IntensityRangeStandardization()
        trained_model, transformed_images = irs.train_transform([i[m] for i, m in zip(images, masks)], surpress_mapping_check = args.ignore)
        logger.info('Saving the trained model as {}...'.format(args.smodel))
        with open(args.smodel, 'wb') as f:
                pickle.dump(trained_model, f)
                
    # save the transformed images
    if args.simages:
        logger.info('Saving intensity transformed images to {}...'.format(args.simages))
        for ti, i, m, h, image_name in zip(transformed_images, images, masks, headers, args.images):
            i[m] = ti
            save(i, '{}/{}'.format(args.simages, image_name.split('/')[-1]), h, args.force)
    
    logger.info('Terminated.')

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)
    
    # check if output image exists (will also be performed before saving, but as the smoothing might be very time intensity, a initial check can save frustration)
    if not args.force:
        if os.path.exists(args.output):
            raise parser.error('The output image {} already exists.'.format(args.output))
    
    # loading image
    data_input, header_input = load(args.input)
    
    # apply the watershed
    logger.info('Applying anisotropic diffusion with settings: niter={} / kappa={} / gamma={}...'.format(args.iterations, args.kappa, args.gamma))
    data_output = anisotropic_diffusion(data_input, args.iterations, args.kappa, args.gamma, get_pixel_spacing(header_input))

    # save file
    save(data_output, args.output, header_input, args.force)
    
    logger.info('Successfully terminated.')