Python pylab.axis() Examples

def plot_confusion_matrix(y_true, y_pred, size=None, normalize=False):
    """plot_confusion_matrix."""
    cm = confusion_matrix(y_true, y_pred)
    fmt = "%d"
    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        fmt = "%.2f"
    xticklabels = list(sorted(set(y_pred)))
    yticklabels = list(sorted(set(y_true)))
    if size is not None:
        plt.figure(figsize=(size, size))
    heatmap(cm, xlabel='Predicted label', ylabel='True label',
            xticklabels=xticklabels, yticklabels=yticklabels,
            cmap=plt.cm.Blues, fmt=fmt)
    if normalize:
        plt.title("Confusion matrix (norm.)")
    else:
        plt.title("Confusion matrix")
    plt.gca().invert_yaxis() 

def plot_roc_curve(y_true, y_score, size=None):
    """plot_roc_curve."""
    false_positive_rate, true_positive_rate, thresholds = roc_curve(
        y_true, y_score)
    if size is not None:
        plt.figure(figsize=(size, size))
        plt.axis('equal')
    plt.plot(false_positive_rate, true_positive_rate, lw=2, color='navy')
    plt.plot([0, 1], [0, 1], color='gray', lw=1, linestyle='--')
    plt.xlabel('False positive rate')
    plt.ylabel('True positive rate')
    plt.ylim([-0.05, 1.05])
    plt.xlim([-0.05, 1.05])
    plt.grid()
    plt.title('Receiver operating characteristic AUC={0:0.2f}'.format(
        roc_auc_score(y_true, y_score))) 

def plot_rectified(self):
        import pylab
        pylab.title('rectified')
        pylab.imshow(self.rectified)

        for line in self.vlines:
            p0, p1 = line
            p0 = self.inv_transform(p0)
            p1 = self.inv_transform(p1)
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='green')

        for line in self.hlines:
            p0, p1 = line
            p0 = self.inv_transform(p0)
            p1 = self.inv_transform(p1)
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='red')

        pylab.axis('image');
        pylab.grid(c='yellow', lw=1)
        pylab.plt.yticks(np.arange(0, self.l, 100.0));
        pylab.xlim(0, self.w)
        pylab.ylim(self.l, 0) 

def click_event(self, event):
		"""Whenever a click occurs on the coefficent axes we modify the coefficents and update the 
plot"""
		
		if event.xdata is None:#we clicked outside the axis
			return

		idx = M.getp(event.inaxes,'label')
		
		print idx, event.xdata, event.ydata
				
		self.c[idx] = event.xdata
		self.c[idx+1] = event.ydata

		self.replotf()
		self.replotc() 

def __init__(self, norder = 2):
		"""Initializes the class when returning an instance. Pass it the polynomial order. It will 
set up two figure windows, one for the graph the other for the coefficent interface. It will then initialize 
the coefficients to zero and plot the (not so interesting) polynomial."""
		
		self.order = norder
		
		self.c = M.zeros(self.order,'f')
		self.ax = [None]*(self.order-1)#M.zeros(self.order-1,'i') #Coefficent axes
		
		self.ffig = M.figure() #The first figure window has the plot
		self.replotf()
		
		self.cfig = M.figure() #The second figure window has the 
		row = M.ceil(M.sqrt(self.order-1))
		for n in xrange(self.order-1):
			self.ax[n] = M.subplot(row, row, n+1)
			M.setp(self.ax[n],'label', n)
			M.plot([0],[0],'.')
			M.axis([-1, 1, -1, 1]);
			
		self.replotc()
		M.connect('button_press_event', self.click_event) 

def transform_to_2d(data, max_axis):
    """
    Projects 3d data cube along one axis using maximum intensity with
    preservation of the signs. Adapted from nilearn.
    """
    import numpy as np

    # get the shape of the array we are projecting to
    new_shape = list(data.shape)
    del new_shape[max_axis]

    # generate a 3D indexing array that points to max abs value in the
    # current projection
    a1, a2 = np.indices(new_shape)
    inds = [a1, a2]
    inds.insert(max_axis, np.abs(data).argmax(axis=max_axis))

    # take the values where the absolute value of the projection
    # is the highest
    maximum_intensity_data = data[inds]

    return np.rot90(maximum_intensity_data) 

def plot_iris_knn():
    iris = datasets.load_iris()
    X = iris.data[:, :2]  # we only take the first two features. We could
                        # avoid this ugly slicing by using a two-dim dataset
    y = iris.target

    knn = neighbors.KNeighborsClassifier(n_neighbors=3)
    knn.fit(X, y)

    x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
    y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.figure()
    pl.pcolormesh(xx, yy, Z, cmap=cmap_light)

    # Plot also the training points
    pl.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    pl.xlabel('sepal length (cm)')
    pl.ylabel('sepal width (cm)')
    pl.axis('tight') 

def plot_iris_knn():
    iris = datasets.load_iris()
    X = iris.data[:, :2]  # we only take the first two features. We could
                        # avoid this ugly slicing by using a two-dim dataset
    y = iris.target

    knn = neighbors.KNeighborsClassifier(n_neighbors=3)
    knn.fit(X, y)

    x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
    y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.figure()
    pl.pcolormesh(xx, yy, Z, cmap=cmap_light)

    # Plot also the training points
    pl.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    pl.xlabel('sepal length (cm)')
    pl.ylabel('sepal width (cm)')
    pl.axis('tight') 

def plot_iris_knn():
    iris = datasets.load_iris()
    X = iris.data[:, :2]  # we only take the first two features. We could
                        # avoid this ugly slicing by using a two-dim dataset
    y = iris.target

    knn = neighbors.KNeighborsClassifier(n_neighbors=3)
    knn.fit(X, y)

    x_min, x_max = X[:, 0].min() - .1, X[:, 0].max() + .1
    y_min, y_max = X[:, 1].min() - .1, X[:, 1].max() + .1
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
                         np.linspace(y_min, y_max, 100))
    Z = knn.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    pl.figure()
    pl.pcolormesh(xx, yy, Z, cmap=cmap_light)

    # Plot also the training points
    pl.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)
    pl.xlabel('sepal length (cm)')
    pl.ylabel('sepal width (cm)')
    pl.axis('tight') 

def plot_original(self):
        import pylab
        pylab.title('original')
        pylab.imshow(self.data)

        for line in self.lines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='blue', alpha=0.3)

        for line in self.vlines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='green')

        for line in self.hlines:
            p0, p1 = line
            pylab.plot((p0[0], p1[0]), (p0[1], p1[1]), c='red')

        pylab.axis('image');
        pylab.grid(c='yellow', lw=1)
        pylab.plt.yticks(np.arange(0, self.l, 100.0));
        pylab.xlim(0, self.w)
        pylab.ylim(self.l, 0) 

def ransac_guess_color(colors, n_iter=50, std=2):
    colors = rgb2lab(colors)
    colors = colors.reshape(-1, 3)
    masked = colors[:, 0] < 0.1
    colors = colors[~masked]
    assert len(colors) > 0, "Must have at least one color"

    best_mu = np.array([0, 0, 0])
    best_n = 0
    for k in range(n_iter):
        subset = colors[np.random.choice(np.arange(len(colors)), 1)]

        mu = subset.mean(0)
        #inliers = (((colors - mu) ** 2 / std) < 1).all(1)
        inliers = ((np.sqrt(np.sum((colors - mu)**2, axis=1))  / std) < 1)

        mu = colors[inliers].mean(0)
        n = len(colors[inliers])
        if n > best_n:
            best_n = n
            best_mu = mu
    #import ipdb; ipdb.set_trace()
    best_mu = np.squeeze(lab2rgb(np.array([[best_mu]])))
    return best_mu 

def padRightDownCorner(img, stride, padValue):
    h = img.shape[0]
    w = img.shape[1]

    pad = 4 * [None]
    pad[0] = 0 # up
    pad[1] = 0 # left
    pad[2] = 0 if (h%stride==0) else stride - (h % stride) # down
    pad[3] = 0 if (w%stride==0) else stride - (w % stride) # right

    img_padded = img
    pad_up = np.tile(img_padded[0:1,:,:]*0 + padValue, (pad[0], 1, 1))
    img_padded = np.concatenate((pad_up, img_padded), axis=0)
    pad_left = np.tile(img_padded[:,0:1,:]*0 + padValue, (1, pad[1], 1))
    img_padded = np.concatenate((pad_left, img_padded), axis=1)
    pad_down = np.tile(img_padded[-2:-1,:,:]*0 + padValue, (pad[2], 1, 1))
    img_padded = np.concatenate((img_padded, pad_down), axis=0)
    pad_right = np.tile(img_padded[:,-2:-1,:]*0 + padValue, (1, pad[3], 1))
    img_padded = np.concatenate((img_padded, pad_right), axis=1)

    return img_padded, pad 

def _cut_windows_vertically(self, door_top, roof_top, sky_sig, win_strip):
        win_sig = np.percentile(win_strip, 85, axis=1)
        win_sig[sky_sig > 0.5] = 0
        if win_sig.max() > 0:
            win_sig /= win_sig.max()
        win_sig[:roof_top] = 0
        win_sig[door_top:] = 0
        runs, starts, values = run_length_encode(win_sig > 0.5)
        win_heights = runs[values]
        win_tops = starts[values]
        if len(win_heights) > 0:
            win_bottom = win_tops[-1] + win_heights[-1]
            win_top = win_tops[0]
            win_vertical_spacing = np.diff(win_tops).mean() if len(win_tops) > 1 else 0
        else:
            win_bottom = win_top = win_vertical_spacing = -1

        self.top = int(win_top)
        self.bottom = int(win_bottom)
        self.vertical_spacing = int(win_vertical_spacing)
        self.vertical_scores = make_list(win_sig)
        self.heights = np.array(win_heights)
        self.tops = np.array(win_tops) 

def spectrum_subplot (spectrum):
    '''Plot a spectrum, with x-axis the wavelength, and y-axis the intensity.
    The curve is colored at that wavelength by the (approximate) color of a
    pure spectral color at that wavelength, with intensity constant over wavelength.
    (This means that dark looking colors here mean that wavelength is poorly viewed by the eye.

    This is not a complete plotting function, e.g. no file is saved, etc.
    It is assumed that this function is being called by one that handles those things.'''
    (num_wl, num_cols) = spectrum.shape
    # get rgb colors for each wavelength
    rgb_colors = numpy.empty ((num_wl, 3))
    for i in range (0, num_wl):
        wl_nm = spectrum [i][0]
        xyz = ciexyz.xyz_from_wavelength (wl_nm)
        rgb_colors [i] = colormodels.rgb_from_xyz (xyz)
    # scale to make brightest rgb value = 1.0
    rgb_max = numpy.max (rgb_colors)
    scaling = 1.0 / rgb_max
    rgb_colors *= scaling
    # draw color patches (thin vertical lines matching the spectrum curve) in color
    for i in range (0, num_wl-1):    # skipping the last one here to stay in range
        x0 = spectrum [i][0]
        x1 = spectrum [i+1][0]
        y0 = spectrum [i][1]
        y1 = spectrum [i+1][1]
        poly_x = [x0,  x1,  x1, x0]
        poly_y = [0.0, 0.0, y1, y0]
        color_string = colormodels.irgb_string_from_rgb (rgb_colors [i])
        pylab.fill (poly_x, poly_y, color_string, edgecolor=color_string)
    # plot intensity as a curve
    pylab.plot (
        spectrum [:,0], spectrum [:,1],
        color='k', linewidth=2.0, antialiased=True) 

def melgram_v2(audio_file_path):
    # Load sound file
    y, sr = librosa.load(audio_file_path)

    # Let's make and display a mel-scaled power (energy-squared) spectrogram
    S = librosa.feature.melspectrogram(y, sr=sr, n_mels=128)

    # Convert to log scale (dB). We'll use the peak power as reference.
    log_S = librosa.core.amplitude_to_db(S, ref=np.max)

    # Make a new figure
    plt.figure(figsize=(12, 4))

    # Display the spectrogram on a mel scale
    # sample rate and hop length parameters are used to render the time axis
    librosa.display.specshow(log_S, sr=sr, x_axis='time', y_axis='mel')

    # Put a descriptive title on the plot
    plt.title('mel power spectrogram')

    # draw a color bar
    plt.colorbar(format='%+02.0f dB')

    # Make the figure layout compact
    plt.tight_layout()
    plt.show() 

def drawModel(mfeat, mode="black", parts=True):
    """
        draw the HOG weight of an object model
    """
    col = ["r", "g", "b"]
    import drawHOG
    lev = len(mfeat)
    if mfeat[0].shape[0] > mfeat[0].shape[1]:
        sy = 1
        sx = lev
    else:
        sy = lev
        sx = 1
    for l in range(lev):
        pylab.subplot(sy, sx, l + 1)
        if mode == "white":
            drawHOG9(mfeat[l])
        elif mode == "black":
            img = drawHOG.drawHOG(mfeat[l])
            pylab.axis("off")
            pylab.imshow(img, cmap=pylab.cm.gray, interpolation="nearest")
        if parts == True:
            for x in range(0, 2 ** l):
                for y in range(0, 2 ** l):
                    boxHOG(mfeat[0].shape[1] * x, mfeat[0].shape[0] * y,
                           mfeat[0].shape[1], mfeat[0].shape[0], col[l], 5 - l) 

def melgram_v1(audio_file_path, to_file):
    sig, fs = librosa.load(audio_file_path)

    pylab.axis('off')  # no axis
    pylab.axes([0., 0., 1., 1.], frameon=False, xticks=[], yticks=[])  # Remove the white edge
    S = librosa.feature.melspectrogram(y=sig, sr=fs)
    librosa.display.specshow(librosa.power_to_db(S, ref=np.max))
    pylab.savefig(to_file, bbox_inches=None, pad_inches=0)
    pylab.close() 

def boxHOG(px, py, dx, dy, col, lw):
    """
        bbox one the HOG weights
    """
    k = 1
    d = 15
    pylab.plot([px * d + 0 - k, px * d + 0 - k],
               [py * d + 0 - k, py * d + dy * d - k], col, lw=lw)
    pylab.plot([px * d + 0 - k, px * d + dx * d - k],
               [py * d + 0 - k, py * d + 0 - k], col, lw=lw)
    pylab.plot([px * d + dx * 15 - k, px * d + dx * d - k],
               [py * d + 0 - k, py * d + dy * d - k], col, lw=lw)
    pylab.plot([px * d + 0 - k, px * d + dx * d - k],
               [py * d + dy * d - k, py * d + dy * d - k], col, lw=lw)
    pylab.axis("image") 

def __list_markings__(self,classification):
        marks_list = self.__classification_to_markings__(classification)
        roi = self.__load_roi__(classification)

        for mark in marks_list:
            x = float(mark["x"])*self.scale
            y = float(mark["y"])*self.scale

            if not("animal" in mark):
                animal_type = None
            else:
                animal_type = mark["animal"]

            #find which line segment on the roi the point lies on (x-axis wise)
            for segment_index in range(len(roi)-1):
                if (roi[segment_index][0] <= x) and (roi[segment_index+1][0] >= x):
                    rX1,rY1 = roi[segment_index]
                    rX2,rY2 = roi[segment_index+1]

                    m = (rY2-rY1)/float(rX2-rX1)
                    rY = m*(x-rX1)+rY1

                    if y >= rY:
                        # we have found a valid marking
                        # create a special type of animal None that is used when the animal type is missing
                        # thus, the marking will count towards not being noise but will not be used when determining the type

                        yield (x,y),animal_type
                        break
                    else:
                        break 

def save_digit(z, filename, cmap=pylab.cm.gray):
    pylab.clf()
    pylab.axis('off')
    pylab.imshow(z.reshape((28, 28)), interpolation='nearest', cmap=cmap, vmin=-1, vmax=1)
    pylab.savefig('results/' + filename + '.pdf')
    pylab.clf() 

def tighten_x_axis (x_list):
    '''Tighten the x axis (only) of the current plot to match the given range of x values.
    The y axis limits are not affected.'''
    x_min = min (x_list)
    x_max = max (x_list)
    pylab.xlim ((x_min, x_max))

#
# Patch plots - Plots with each color value as a solid patch, with optional labels.
# 

def plot_polynomial_regression():
    rng = np.random.RandomState(0)
    x = 2*rng.rand(100) - 1

    f = lambda t: 1.2 * t**2 + .1 * t**3 - .4 * t **5 - .5 * t ** 9
    y = f(x) + .4 * rng.normal(size=100)

    x_test = np.linspace(-1, 1, 100)

    pl.figure()
    pl.scatter(x, y, s=4)

    X = np.array([x**i for i in range(5)]).T
    X_test = np.array([x_test**i for i in range(5)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='4th order')

    X = np.array([x**i for i in range(10)]).T
    X_test = np.array([x_test**i for i in range(10)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='9th order')

    pl.legend(loc='best')
    pl.axis('tight')
    pl.title('Fitting a 4th and a 9th order polynomial')

    pl.figure()
    pl.scatter(x, y, s=4)
    pl.plot(x_test, f(x_test), label="truth")
    pl.axis('tight')
    pl.title('Ground truth (9th order polynomial)') 

def rgb_patch_plot (
    rgb_colors,
    color_names,
    title,
    filename,
    patch_gap = 0.05,
    num_across = 6):
    '''Draw a set of color patches, specified as linear rgb colors.'''

    def draw_patch (x0, y0, color, name, patch_gap):
        '''Draw a patch of color.'''
        # patch relative vertices
        m = patch_gap
        omm = 1.0 - m
        poly_dx = [m, m, omm, omm]
        poly_dy = [m, omm, omm, m]
        # construct vertices
        poly_x = [ x0 + dx_i for dx_i in poly_dx ]
        poly_y = [ y0 + dy_i for dy_i in poly_dy ]
        pylab.fill (poly_x, poly_y, color)
        if name != None:
            dtext = 0.1
            pylab.text (x0+dtext, y0+dtext, name, size=8.0)

    # make plot with each color with one patch
    pylab.clf()
    num_colors = len (rgb_colors)
    for i in range (0, num_colors):
        (iy, ix) = divmod (i, num_across)
        # get color as a displayable string
        colorstring = colormodels.irgb_string_from_rgb (rgb_colors [i])
        if color_names != None:
            name = color_names [i]
        else:
            name = None
        draw_patch (float (ix), float (-iy), colorstring, name, patch_gap)
    pylab.axis ('off')
    pylab.title (title)
    print ('Saving plot %s' % str (filename))
    pylab.savefig (filename) 

def plot_polynomial_regression():
    rng = np.random.RandomState(0)
    x = 2*rng.rand(100) - 1

    f = lambda t: 1.2 * t**2 + .1 * t**3 - .4 * t **5 - .5 * t ** 9
    y = f(x) + .4 * rng.normal(size=100)

    x_test = np.linspace(-1, 1, 100)

    pl.figure()
    pl.scatter(x, y, s=4)

    X = np.array([x**i for i in range(5)]).T
    X_test = np.array([x_test**i for i in range(5)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='4th order')

    X = np.array([x**i for i in range(10)]).T
    X_test = np.array([x_test**i for i in range(10)]).T
    regr = linear_model.LinearRegression()
    regr.fit(X, y)
    pl.plot(x_test, regr.predict(X_test), label='9th order')

    pl.legend(loc='best')
    pl.axis('tight')
    pl.title('Fitting a 4th and a 9th order polynomial')

    pl.figure()
    pl.scatter(x, y, s=4)
    pl.plot(x_test, f(x_test), label="truth")
    pl.axis('tight')
    pl.title('Ground truth (9th order polynomial)') 

def replotc(self):
		"""This replots the coefficients."""
	
		M.figure(self.cfig.number)
		M.ioff()
		for n in xrange(self.order-1):
			M.axes(self.ax[n])
			#M.cla()
			M.plot([self.c[n]], [self.c[n+1]],'ko')
			M.xlabel("$c_%d$" %(n))
			M.ylabel("$c_%d$" %(n+1))
			M.axis([-1, 1, -1, 1]);
			del self.ax[n].lines[0]

		M.draw() 

def melgram_v1(audio_file_path, to_file):
    sig, fs = librosa.load(audio_file_path)

    pylab.axis('off')  # no axis
    pylab.axes([0., 0., 1., 1.], frameon=False, xticks=[], yticks=[])  # Remove the white edge
    S = librosa.feature.melspectrogram(y=sig, sr=fs)
    librosa.display.specshow(librosa.power_to_db(S, ref=np.max))
    pylab.savefig(to_file, bbox_inches=None, pad_inches=0)
    pylab.close() 

def melgram_v2(audio_file_path):
    # Load sound file
    y, sr = librosa.load(audio_file_path)

    # Let's make and display a mel-scaled power (energy-squared) spectrogram
    S = librosa.feature.melspectrogram(y, sr=sr, n_mels=128)

    # Convert to log scale (dB). We'll use the peak power as reference.
    log_S = librosa.core.amplitude_to_db(S, ref=np.max)

    # Make a new figure
    plt.figure(figsize=(12, 4))

    # Display the spectrogram on a mel scale
    # sample rate and hop length parameters are used to render the time axis
    librosa.display.specshow(log_S, sr=sr, x_axis='time', y_axis='mel')

    # Put a descriptive title on the plot
    plt.title('mel power spectrogram')

    # draw a color bar
    plt.colorbar(format='%+02.0f dB')

    # Make the figure layout compact
    plt.tight_layout()
    plt.show() 

def create_figure():
    psd = test_pcovar()
    pylab.axis([-0.5,0.5,-60,0])
    pylab.savefig('psd_covar.png') 

def plot_facade_cuts(self):

        facade_sig = self.facade_edge_scores.sum(0)
        facade_cuts = find_facade_cuts(facade_sig, dilation_amount=self.facade_merge_amount)
        mu = np.mean(facade_sig)
        sigma = np.std(facade_sig)

        w = self.rectified.shape[1]
        pad=10

        gs1 = pl.GridSpec(5, 5)
        gs1.update(wspace=0.5, hspace=0.0)  # set the spacing between axes.

        pl.subplot(gs1[:3, :])
        pl.imshow(self.rectified)
        pl.vlines(facade_cuts, *pl.ylim(), lw=2, color='black')
        pl.axis('off')
        pl.xlim(-pad, w+pad)

        pl.subplot(gs1[3:, :], sharex=pl.gca())
        pl.fill_between(np.arange(w), 0, facade_sig, lw=0, color='red')
        pl.fill_between(np.arange(w), 0, np.clip(facade_sig, 0, mu+sigma), color='blue')
        pl.plot(np.arange(w), facade_sig, color='blue')

        pl.vlines(facade_cuts, facade_sig[facade_cuts], pl.xlim()[1], lw=2, color='black')
        pl.scatter(facade_cuts, facade_sig[facade_cuts])

        pl.axis('off')

        pl.hlines(mu, 0, w, linestyle='dashed', color='black')
        pl.text(0, mu, '$\mu$ ', ha='right')

        pl.hlines(mu + sigma, 0, w, linestyle='dashed', color='gray',)
        pl.text(0, mu + sigma, '$\mu+\sigma$ ', ha='right')
        pl.xlim(-pad, w+pad) 

def _solve_lr(vlines, w, l, opt_options=OPTIMIZATION_OPTIONS, opt_method=OPTIMIZATION_METHOD, limit=0.3):
    """ Solve for the left and right edge displacement.
    This routine estimates the amount to move the upper left and right cornders of the image
     in a horizontal direction in order to make the given lines parallel and vertical.

    :param vlines: Lines that we want to map to vertical lines.
    :param w:  The width of the image
    :param l:   The height of the image
    :param opt_options: Optimization options passed into `minimize`
    :param opt_method: The optimization method.
    :param limit:  A limit on the amount of displacement -- beyond this and we will assume failure.
    :return: (dl, dr),   the horizontal displacement of the left and right corners.
    """
    if len(vlines) == 0:
        return 0, 0

    a = np.append(vlines[:, 0, :], np.ones((len(vlines), 1)), axis=1)
    b = np.append(vlines[:, 1, :], np.ones((len(vlines), 1)), axis=1)

    def objective(x):
        dl, dr = x
        Hv = np.linalg.inv(H_v(dl, dr, w, l))
        return np.sum(np.abs(Hv[0, :].dot(a.T) / Hv[2, :].dot(a.T) - Hv[0, :].dot(b.T) / Hv[2, :].dot(b.T)))

    res = minimize(objective, (0., 0.),
                   options=opt_options,
                   method=opt_method)
    dl, dr = res.x

    # Give up if the solution is not plausible (this indicates that the 'vlines' are too noisy
    if abs(dl) > limit * w:
        dl = 0
    if abs(dr) > limit * w:
        dr = 0
    return dl, dr