def METHOD_NAME():
    parser = _build_arg_parser()
    args = parser.parse_args()
    if args.verbose:
        logging.getLogger().setLevel(logging.INFO)
    inputs = [args.in_bvec, args.in_peaks, args.in_FA]
    optional = [args.mask, args.peaks_vals]
    assert_inputs_exist(parser, inputs, optional=optional)
    assert_outputs_exist(parser, args, args.out_bvec)
    _, bvecs = read_bvals_bvecs(None, args.in_bvec)
    fa = nib.load(args.in_FA).get_fdata()
    peaks = nib.load(args.in_peaks).get_fdata()
    # convert peaks to a volume of shape (H, W, D, N, 3)
    if args.column_wise:
        peaks = np.reshape(peaks, peaks.shape[:3] + (3, -1))
        peaks = np.transpose(peaks, axes=(0, 1, 2, 4, 3))
    else:
        peaks = np.reshape(peaks, peaks.shape[:3] + (-1, 3))
    N = peaks.shape[3]
    if N > 1:
        if not args.peaks_vals:
            parser.error('More than one principal direction per voxel. Specify'
                         ' peaks values with --peaks_vals to proceed.')
        peaks_vals = nib.load(args.peaks_vals).get_fdata()
        indices_max = np.argmax(peaks_vals, axis=-1)[..., None, None]
        peaks = np.take_along_axis(peaks, indices_max, axis=-2)
    peaks = np.squeeze(peaks)
    if args.mask:
        mask = get_data_as_mask(nib.load(args.mask))
        fa[np.logical_not(mask)] = 0
        peaks[np.logical_not(mask)] = 0
    peaks[fa < args.fa_th] = 0
    coherence, transform =\
        compute_fiber_coherence_table(peaks, fa)
    best_t = transform[np.argmax(coherence)]
    if (best_t == np.eye(3)).all():
        logging.info('b-vectors are already correct.')
        correct_bvecs = bvecs
    else:
        logging.info('Applying correction to b-vectors. '
                     'Transform is: \n{0}.'.format(best_t))
        correct_bvecs = np.dot(bvecs, best_t)
    logging.info('Saving bvecs to file: {0}.'.format(args.out_bvec))
    # FSL format (3, N)
    np.savetxt(args.out_bvec, correct_bvecs.T, '%.8f')