def calc_hist(hsv,mask, figure=None, ax1=None, ax2=None):
    chans = cv2.split(hsv)
    # Or maybe I should use Numpy indexing for faster splitting: h=hsv[:,:,0] 
    hist_h = cv2.calcHist([chans[0]], [0], mask, [180], [0, 180])
    hist_s = cv2.calcHist([chans[1]], [0], mask, [256], [0, 256])
    #print hist_s
    hist_h = hist_h.flatten()
    hist_s = hist_s.flatten()
    # Apply Gaussian low pass for histogram (using scipy)
    hist_hg = ndi.gaussian_filter1d(hist_h, sigma=1.5, output=np.float64, mode='nearest')
    hist_sg = ndi.gaussian_filter1d(hist_s, sigma=1.5, output=np.float64, mode='nearest')
    hue_max = np.argmax(hist_hg)
    saturation_max = np.argmax(hist_sg) if np.argmax(hist_sg) >= 20 else 20
    #print hue_max, saturation_max
    #ax1.clear(), ax2.clear()
    #ax1.set_autoscale_on(False)
    # ax1.plot(range(180),hist_hg)
    # ax2.plot(range(256),hist_sg)
    # ax1.set_ylim([0,1200])
    # ax1.set_xlim([0,180])
    # ax2.set_xlim([0,256])
    # ax2.set_ylim([0,1200])
    # figure.canvas.draw()
    #plt.xlim([0, 180])
    lower = np.array([hue_max+20,saturation_max-20,20])
    upper = np.array([hue_max+20,saturation_max+20,255])
    mask_color = cv2.inRange(hsv, lower, upper)
    return hue_max, hist_hg, saturation_max, hist_sg, mask_color

def _smooth_array(arr, affine, fwhm=None, ensure_finite=True, copy=True):
    """Smooth images by applying a Gaussian filter.

    Apply a Gaussian filter along the three first dimensions of arr.

    Parameters
    ==========
    arr: numpy.ndarray
        4D array, with image number as last dimension. 3D arrays are also
        accepted.

    affine: numpy.ndarray
        (4, 4) matrix, giving affine transformation for image. (3, 3) matrices
        are also accepted (only these coefficients are used).

    fwhm: scalar or numpy.ndarray
        Smoothing strength, as a full-width at half maximum, in millimeters.
        If a scalar is given, width is identical on all three directions.
        A numpy.ndarray must have 3 elements, giving the FWHM along each axis.
        If fwhm is None, no filtering is performed (useful when just removal
        of non-finite values is needed)

    ensure_finite: bool
        if True, replace every non-finite values (like NaNs) by zero before
        filtering.

    copy: bool
        if True, input array is not modified. False by default: the filtering
        is performed in-place.

    Returns
    =======
    filtered_arr: numpy.ndarray
        arr, filtered.

    Notes
    =====
    This function is most efficient with arr in C order.
    """

    if copy:
        arr = arr.copy()

    # Keep only the scale part.
    affine = affine[:3, :3]

    if ensure_finite:
        # SPM tends to put NaNs in the data outside the brain
        arr[np.logical_not(np.isfinite(arr))] = 0

    if fwhm is not None:
        # Convert from a FWHM to a sigma:
        # Do not use /=, fwhm may be a numpy scalar
        fwhm = fwhm / np.sqrt(8 * np.log(2))
        vox_size = np.sqrt(np.sum(affine ** 2, axis=0))
        sigma = fwhm / vox_size
        for n, s in enumerate(sigma):
            ndimage.gaussian_filter1d(arr, s, output=arr, axis=n)

    return arr

    def _homogenize_params(self, other, maxdiff=1):
        """
        Return triple with a tuple of indices (in self and other, respectively),
        factors and constants at these frequencies.
        
        Parameters
        ----------
        other : CallistoSpectrogram
            Spectrogram to be homogenized with the current one.
        maxdiff : float
            Threshold for which frequencies are considered equal.
        """

        pairs_indices = [(x, y) for x, y, d in minimal_pairs(self.freq_axis, other.freq_axis) if d <= maxdiff]

        pairs_data = [(self[n_one, :], other[n_two, :]) for n_one, n_two in pairs_indices]

        # XXX: Maybe unnecessary.
        pairs_data_gaussian = [(gaussian_filter1d(a, 15), gaussian_filter1d(b, 15)) for a, b in pairs_data]

        # If we used integer arithmetic, we would accept more invalid
        # values.
        pairs_data_gaussian64 = np.float64(pairs_data_gaussian)
        least = [leastsq(self._to_minimize(a, b), [1, 0])[0] for a, b in pairs_data_gaussian64]

        factors = [x for x, y in least]
        constants = [y for x, y in least]

        return pairs_indices, factors, constants

def test_multiple_modes_sequentially():
    # Test that the filters with multiple mode cababilities for different
    # dimensions give the same result as applying the filters with
    # different modes sequentially
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    modes = ['reflect', 'wrap']

    expected = sndi.gaussian_filter1d(arr, 1, axis=0, mode=modes[0])
    expected = sndi.gaussian_filter1d(expected, 1, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.gaussian_filter(arr, 1, mode=modes))

    expected = sndi.uniform_filter1d(arr, 5, axis=0, mode=modes[0])
    expected = sndi.uniform_filter1d(expected, 5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.uniform_filter(arr, 5, mode=modes))

    expected = sndi.maximum_filter1d(arr, size=5, axis=0, mode=modes[0])
    expected = sndi.maximum_filter1d(expected, size=5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.maximum_filter(arr, size=5, mode=modes))

    expected = sndi.minimum_filter1d(arr, size=5, axis=0, mode=modes[0])
    expected = sndi.minimum_filter1d(expected, size=5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.minimum_filter(arr, size=5, mode=modes))

def bootdensity(data, min, max, nboot, ci):
    """ Calculate density and confidence intervals on density
    for a 1D array of points.  Bandwidth is selected automatically.
    """
    r("""
      limdensity <- function(data, weights=NULL, bw="nrd0")
      {
        density(data, from=%f, to=%f, weights=weights, bw=bw)
      }
      """%(min, max))
    density = r.limdensity(data)
    xdens = N.array(density['x'])
    ydens = N.array(density['y'])
    bw = density['bw']
    #print 'bandwidth:', bw
    ydensboot = N.zeros((nboot, len(xdens)), N.float)
    ndata = len(data)
    ran = N.random.uniform(0, ndata, (nboot,ndata)).astype(N.int)
    for i in range(nboot):
        den = r.limdensity(data[ran[i]])
        y = N.array(den['y'])
        ydensboot[i] = y
    ydensbootsort = N.sort(ydensboot, axis=0)
    ydensbootsort = interp1d(N.arange(0, 1.000001, 1.0/(nboot-1)),
                             ydensbootsort, axis=0)
    ilow = (0.5-ci/2.0)
    ihigh = (0.5+ci/2.0)
    ydenslow, ydenshigh = ydensbootsort((ilow, ihigh))
    ydenslow = gaussian_filter1d(ydenslow, bw*512/10.0)
    ydenshigh, ydenshigh = ydensbootsort((ihigh, ihigh))
    ydenshigh = gaussian_filter1d(ydenshigh, bw*512/10.0)
    return xdens, ydens, ydenslow, ydenshigh, bw

def createOrdered3D(centerPoints,edgeSets):
	f = open(name + 'TrackPointsTest.xyz', 'w')
	fig = plt.figure()
	ax = fig.add_subplot(111, projection='3d')
	ax.set_zbound(lower=0, upper=1400)
	for edgeSet in edgeSets:
		xArr = []
		yArr = []
		zArr = []
		for point in edgeSet:
			xArr.append(point[0])
			yArr.append(point[1])
			height = getHeight(point[0],point[1],centerPoints)*HEIGHT_SCALE
			zArr.append(height)
			f.write(str(point[0]) + " "+ str(point[1]) + " " + str(height)+"\n")
		t = np.linspace(0, 1, len(zArr))
		t2 = np.linspace(0, 1, len(zArr))

		x2 = np.interp(t2, t, xArr)
		y2 = np.interp(t2, t, yArr)
		z2 = np.interp(t2, t, zArr)
		sigma = 10
		x3 = gaussian_filter1d(x2, sigma)
		y3 = gaussian_filter1d(y2, sigma)
		z3 = gaussian_filter1d(z2, sigma)
		ax.plot(xArr,yArr,zArr,color = 'b')
	plt.show()

def test_orders_gauss():
    # Check order inputs to Gaussians
    arr = np.zeros((1,))
    assert_equal(0, sndi.gaussian_filter(arr, 1, order=0))
    assert_equal(0, sndi.gaussian_filter(arr, 1, order=3))
    assert_raises(ValueError, sndi.gaussian_filter, arr, 1, -1)
    assert_equal(0, sndi.gaussian_filter1d(arr, 1, axis=-1, order=0))
    assert_equal(0, sndi.gaussian_filter1d(arr, 1, axis=-1, order=3))
    assert_raises(ValueError, sndi.gaussian_filter1d, arr, 1, -1, -1)

 def extract_arrays(self):
     self.xs, self.ys = self.chromatogram.as_arrays()
     if self.smooth:
         self.ys = gaussian_filter1d(self.ys, 1)
     if len(self.xs) > MAX_POINTS:
         new_xs = np.linspace(self.xs.min(), self.xs.max(), MAX_POINTS)
         new_ys = np.interp(new_xs, self.xs, self.ys)
         self.xs = new_xs
         self.ys = new_ys
         self.ys = gaussian_filter1d(self.ys, 1)

def fit_semiconductor(t, data, sav_n=11, sav_deg=4, mode='sav', tr=0.4):
    from scipy.signal import savgol_filter
    from scipy.ndimage import gaussian_filter1d
    from scipy.optimize import leastsq
    ger =   data[..., -1].sum(2).squeeze()
    plt.subplot(121)
    plt.title('Germanium sum')
    plt.plot(t, ger[:,  0])
    plt.plot(t, ger[:,  1])
    if mode =='sav':
        plt.plot(t, savgol_filter(ger[:, 0], sav_n, sav_deg, 0))
        plt.plot(t, savgol_filter(ger[:, 1], sav_n, sav_deg, 0))
    plt.xlim(-1, 3)
    plt.subplot(122)
    plt.title('First dervitate')
    if mode == 'sav':
        derv0 = savgol_filter(ger[:, 0], sav_n, sav_deg, 1)
        derv1 = savgol_filter(ger[:, 1], sav_n, sav_deg, 1)
    elif mode == 'gauss':
        derv0 =  gaussian_filter1d(ger[:, 0], sav_n, order=1)
        derv1 =  gaussian_filter1d(ger[:, 1], sav_n, order=1)
    plt.plot(t , derv0)
    plt.plot(t , derv1)
    plt.xlim(-.8, .8)
    plt.ylim(0, 700)
    plt.minorticks_on()
    plt.grid(1)

    def gaussian(p, ch, res=True):

        i, j = dv.fi(t, -tr), dv.fi(t, tr)
        w = p[0]
        A = p[1]
        x0 = p[2]
        fit = A*np.exp(-(t[i:j]-x0)**2/(2*w**2))
        if res:
            return fit-ch[i:j]
        else:
            return fit


    x0 = leastsq(gaussian, [.2, max(derv0), 0], derv0)
    plt.plot(t[dv.fi(t, -tr):dv.fi(t, tr)], gaussian(x0[0], 0, 0), '--k', )
    plt.text(0.05, 0.9, 'x$_0$ = %.2f\nFWHM = %.2f\nA = %.1f\n'%(x0[0][2],2.35*x0[0][0], x0[0][1]),
             transform=plt.gca().transAxes, va='top')

    x0 = leastsq(gaussian, [.2, max(derv1), 0], derv1)
    plt.plot(t[dv.fi(t, -tr):dv.fi(t, tr)], gaussian(x0[0], 1, 0), '--b', )

    plt.xlim(-.8, .8)
    plt.minorticks_on()
    plt.grid(0)
    plt.tight_layout()
    plt.text(0.5, 0.9, 'x$_0$ = %.2f\nFWHM = %.2f\nA = %.1f\n'%(x0[0][2],2.35*x0[0][0], x0[0][1]),
             transform=plt.gca().transAxes, va='top')

def trace_profile(image, sigma=5., width_factor=1., check_vertical=False):
    """Trace the intensity profile of a tubular structure in an image.

    Parameters
    ----------
    image : array of int or float, shape (M, N[, P])
        The input image. If 3D, the first dimension is flattened by
        summing along that axis.
    sigma : float, optional
        Convolve the intensity with this sigma to estimate the start
        and end of the scan lines.
    width_factor : float, optional
        The width of the line profile is determined automatically, then
        multiplied by this factor.
    check_vertical : bool, optional
        Check whether the tube is arranged top-to-bottom in the image.
        If `False`, it is assumed to be vertical, otherwise, the
        orientation is automatically determined from the image.

    Returns
    -------
    profile : 1D array of float
        The intensity profile of the tube.

    Examples
    --------
    >>> edges = np.array([8, 16, 22, 16, 8])
    >>> middle = np.array([0, 0, 0, 0, 0])
    >>> image = np.vstack([edges, middle, edges])
    >>> trace_profile(image, sigma=0)
    array([ 18.,   0.,  18.])
    >>> image3d = np.array([image, image, image])
    >>> trace_profile(image3d, sigma=0)
    array([ 54.,   0.,  54.])
    >>> trace_profile(image.T, sigma=0, check_vertical=True)
    array([ 18.,   0.,  18.])
    """
    if image.ndim > 2:
        image = image.sum(axis=0)
    if check_vertical:
        top_bottom_mean = np.mean(image[[0, image.shape[0] - 1], :])
        left_right_mean = np.mean(image[:, [0, image.shape[1] - 1]])
        if top_bottom_mean < left_right_mean:
            image = image.T
    top_distribution = nd.gaussian_filter1d(image[0], sigma)
    bottom_distribution = nd.gaussian_filter1d(image[-1], sigma)
    top_loc, top_whm = estimate_mode_width(top_distribution)
    bottom_loc, bottom_whm = estimate_mode_width(bottom_distribution)
    angle = np.arctan(np.abs(float(bottom_loc - top_loc)) / image.shape[0])
    width = np.int(np.ceil(max(top_whm, bottom_whm) * np.cos(angle)))
    profile = profile_line(image,
                           (0, top_loc), (image.shape[0] - 1, bottom_loc),
                           linewidth=width, mode='nearest')
    return profile

def derivatives(flux, dx, s_factor):
    dxdxdx = dx * dx * dx
    # First derivative
    gf = ndimage.gaussian_filter1d(flux, sigma=s_factor, order=1, mode='wrap') / dx

    # Second derivative
    ggf = ndimage.gaussian_filter1d(flux, sigma=s_factor, order=2, mode='wrap') / (dx * dx)

    # Third derivative
    gggf = np.array(ndimage.gaussian_filter1d(flux, sigma=s_factor, order=3, mode='wrap') / dxdxdx)

    return gf, ggf, gggf

def computeSegmentGradients(segmentArray):
    xs, ys, zs, ii, ij = segmentArray
    
    xd = ndimage.gaussian_filter1d(xs, 1.0, order=1)
    yd = ndimage.gaussian_filter1d(ys, 1.0, order=1)
    zd = ndimage.gaussian_filter1d(zs, 1.0, order=1)
    
    mag = sqrt( xd**2 + yd**2 + zd**2 )
    mag[mag==0] = 1
    

    return c_[xd/mag, yd/mag, zd/mag].T

	def getBoundsInfo(self):
		'''
		Each variable needs an upper and a lower bound. This portion of the
		algorithm needs some fine-tuning... these bounds determine the feasible
		set. Early experiments show that defining a tight feasible set
		produces quite good results, so that is why the multipliers here 
		are quite small.

		The primal variables are ordered [mu, theta], and the constraints are
		order [L, R, S, T]. The fast index is over position (i.e. there are
		nWrap entries before the variables changes).
		'''

		# Extract single columns
		L = self.columns['%04d.target' % (self.column)]['L']
		R = self.columns['%04d.target' % (self.column)]['R']
		S = self.columns['%04d.target' % (self.column)]['S']
		T = self.columns['%04d.target' % (self.column)]['T']

		# Bounds
		lowerBound = np.zeros(self.optDict['nPrimal'], dtype=float)
		lowerBound[self.optDict['nWrap']:] = 0.33
		upperBound = np.empty(self.optDict['nPrimal'], dtype=float)
		upperBound[:self.optDict['nWrap']] = 30.0 * self.parameter_scale
		upperBound[self.optDict['nWrap']:] = 0.34 * self.parameter_scale #3.0 / 4.0

		# Constraints
		factor = 1e-2
		lowerConstraint = np.empty(self.optDict['nConstraint'], dtype=float)
		upperConstraint = np.empty(self.optDict['nConstraint'], dtype=float)
		for i, val in enumerate([L, R, S, T]):
			lowerConstraint[i*self.optDict['nWrap']:(i+1)*self.optDict['nWrap']] = \
			ndimage.gaussian_filter1d((val - factor * val), 10)
		for i, val in enumerate([L, R, S, T]):
			upperConstraint[i*self.optDict['nWrap']:(i+1)*self.optDict['nWrap']] = \
			ndimage.gaussian_filter1d((val + factor * val), 10)
		for i, name in enumerate(['L','R','S','T']):
			self.sliceDict['con%s' % (name)] = slice(
				i*self.optDict['nWrap'],(i+1)*self.optDict['nWrap'])

		# Save these values in a dict
		self.constraintDict['lowerBound'] = lowerBound
		self.constraintDict['upperBound'] = upperBound
		self.constraintDict['lowerConstraint'] = lowerConstraint
		self.constraintDict['upperConstraint'] = upperConstraint
		self.constraintDict['constraintArray'] = np.empty_like(lowerConstraint)

		# Get the model update from FWI
		self.fwiTheta = self.columns['%04d.gradient' % (self.column)]['theta']
		self.fwiMu = self.columns['%04d.gradient' % (self.column)]['mu']

		self.fwiTheta *= self.parameter_scale #ndimage.filters.gaussian_filter1d(self.fwiTheta, 10)

 def induced_voltage_generation(self, Beam, length = 'slice_frame'):
     '''
     *Method to calculate the induced voltage through the derivative of the
     profile; the impedance must be of inductive type.*
     '''
     
     index = self.current_turn[0]
         
     if self.periodicity:
         self.derivative_line_density_not_filtered = np.zeros(self.slices.n_slices)
         find_index_slice = np.searchsorted(self.slices.edges, self.t_rev[index])
         if self.smooth_before_after[0]:
             if self.filter_ind_imp == 'gaussian':
                 self.slices.n_macroparticles = ndimage.gaussian_filter1d(self.slices.n_macroparticles, sigma=self.filter_options, mode='wrap')
             elif self.filter_ind_imp == 'chebyshev':
                 nCoefficients, b, a = self.slices.beam_profile_filter_chebyshev(self.filter_options)
             else:
                 raise RuntimeError('filter method not recognised')
         if (self.t_rev[index]-self.slices.bin_centers[find_index_slice])>0: 
             temp = np.concatenate((np.array([self.slices.n_macroparticles[find_index_slice]]), self.slices.n_macroparticles[:find_index_slice+1], np.array([self.slices.n_macroparticles[0]])))
         else:
             temp = np.concatenate((np.array([self.slices.n_macroparticles[find_index_slice-1]]), self.slices.n_macroparticles[:find_index_slice], self.slices.n_macroparticles[:2]))
         self.derivative_line_density_not_filtered[: find_index_slice+1] = np.gradient(temp, self.slices.bin_centers[1]-self.slices.bin_centers[0])[1:-1] / (self.slices.bin_centers[1] - self.slices.bin_centers[0])
         if self.smooth_before_after[1]:
             if self.filter_ind_imp == 'gaussian':
                 self.derivative_line_density_filtered = ndimage.gaussian_filter1d(self.derivative_line_density_not_filtered, sigma=self.filter_options, mode='wrap')
             elif self.filter_ind_imp == 'chebyshev':
                 self.derivative_line_density_filtered = filtfilt(b, a, self.derivative_line_density_not_filtered)
                 self.derivative_line_density_filtered = np.ascontiguousarray(self.derivative_line_density_filtered)
             else:
                 raise RuntimeError('filter method not recognised')
             induced_voltage = - Beam.charge * e * Beam.ratio * \
             self.Z_over_n[0][index] * \
             self.derivative_line_density_filtered / (2 * np.pi * self.revolution_frequency[index])   
         else:
             induced_voltage = - Beam.charge * e * Beam.ratio * \
             self.Z_over_n[0][index] * \
             self.derivative_line_density_not_filtered / (2 * np.pi * self.revolution_frequency[index])
     else:
         induced_voltage = - Beam.charge * e / (2 * np.pi) * Beam.ratio * \
             self.Z_over_n[0][index] / self.revolution_frequency[index] * \
             self.slices.beam_profile_derivative(self.deriv_mode)[1] / \
             (self.slices.bin_centers[1] - self.slices.bin_centers[0])    
     
     self.induced_voltage = induced_voltage[0:self.slices.n_slices]
     
     if isinstance(length, int):
         max_length = len(induced_voltage)
         if length > max_length:
             induced_voltage = np.lib.pad(self.induced_voltage, (0, length - max_length), 'constant', constant_values=(0,0))
         return induced_voltage[0:length]

def handler(points, mr, gofscale, gof, sigma):
    from pdf2py import readwrite
    from meg import density
    from mri import transform
    from scipy import ndimage
    from nifti import NiftiImage
    from numpy import float32, int16, array

    report = {}
    fids = eval(mr.description)
    lpa = fids[0]
    rpa = fids[1]
    nas = fids[2]
    # self.points = array([[0,0,0],[10,0,0],[0,20,0]])#DEBUG-----------------
    xyz = transform.meg2mri(lpa, rpa, nas, dipole=points)
    # readwrite.writedata(xyz, os.path.dirname(mripath)+'/'+'xyz')
    print "lpa, rpa, nas", lpa, rpa, nas
    print mr.pixdim

    # do some scaling of the dips using the GOF as a weight.
    VoxDim = mr.voxdim[::-1]
    xyzscaled = (xyz / VoxDim).T
    print xyzscaled
    d = density.calc(xyz)
    gofscale = float32(gofscale)
    print "gofscale", gofscale
    s = gof - gofscale
    sf = (1 / (1 - gofscale)) * s
    ds = d * sf

    # apply a 1D gaussian filter
    z = density.val2img(mr.data, ds, xyzscaled)
    # sigma = float32(self.sigmaval.GetValue())
    print "sigma", sigma
    # sigma = 3
    print "filtering 1st dimension"
    f = ndimage.gaussian_filter1d(z, sigma * 1 / VoxDim[0], axis=0)
    print "filtering 2nd dimension"
    f = ndimage.gaussian_filter1d(f, sigma * 1 / VoxDim[1], axis=1)
    print "filtering 3rd dimension"
    f = ndimage.gaussian_filter1d(f, sigma * 1 / VoxDim[2], axis=2)

    scaledf = int16((z.max() / f.max()) * f * 1000)
    print "writing nifti output image"
    overlay = NiftiImage(int16(scaledf))

    overlay.setDescription(mr.description)
    overlay.setFilename(mr.filename + "dd")
    overlay.setQForm(mr.getQForm())

    return overlay

def _smooth_data_array(arr, affine, fwhm, copy=True):
    """Smooth images with a a Gaussian filter.

    Apply a Gaussian filter along the three first dimensions of arr.

    Parameters
    ----------
    arr: numpy.ndarray
        3D or 4D array, with image number as last dimension.

    affine: numpy.ndarray
        Image affine transformation matrix for image.

    fwhm: scalar, numpy.ndarray
        Smoothing kernel size, as Full-Width at Half Maximum (FWHM) in millimeters.
        If a scalar is given, kernel width is identical on all three directions.
        A numpy.ndarray must have 3 elements, giving the FWHM along each axis.

    copy: bool
        if True, will make a copy of the input array. Otherwise will directly smooth the input array.

    Returns
    -------
    smooth_arr: numpy.ndarray
    """

    if arr.dtype.kind == 'i':
        if arr.dtype == np.int64:
            arr = arr.astype(np.float64)
        else:
            arr = arr.astype(np.float32)
    if copy:
        arr = arr.copy()

    # Zeroe possible NaNs and Inf in the image.
    arr[np.logical_not(np.isfinite(arr))] = 0

    try:
        # Keep the 3D part of the affine.
        affine = affine[:3, :3]

        # Convert from FWHM in mm to a sigma.
        fwhm_sigma_ratio = np.sqrt(8 * np.log(2))
        vox_size         = np.sqrt(np.sum(affine ** 2, axis=0))
        sigma            = fwhm / (fwhm_sigma_ratio * vox_size)
        for n, s in enumerate(sigma):
            ndimage.gaussian_filter1d(arr, s, output=arr, axis=n)
    except:
        raise ValueError('Error smoothing the array.')
    else:
        return arr

def gaussian_filter():
    shape = (40, 41, 42, 100)
    smooth_sigma = np.asarray((1., 2, 3.))

    data1 = random_matrix(shape, order='F')
    data2 = order_preserving_copy(data1)

    for img in np.rollaxis(data1, -1):
        img[...] = ndimage.gaussian_filter(img, smooth_sigma)

    for n, s in enumerate(smooth_sigma):
        ndimage.gaussian_filter1d(data2, s, output=data2, axis=n)

    np.testing.assert_almost_equal(data1, data2)

def get_y_derivativemap(flat, flat_bpix, bg_std_norm,
                        max_sep_order=150, pad=50,
                        med_filter_size=(7, 7),
                        flat_mask=None):

    """
    flat
    flat_bpix : bpix'ed flat
    """

    # 1d-derivatives along y-axis : 1st attempt
    # im_deriv = ni.gaussian_filter1d(flat, 1, order=1, axis=0)

    # 1d-derivatives along y-axis : 2nd attempt. Median filter first.

    flat_deriv_bpix = ni.gaussian_filter1d(flat_bpix, 1,
                                           order=1, axis=0)

    # We also make a median-filtered one. This one will be used to make masks.
    flat_medianed = ni.median_filter(flat,
                                     size=med_filter_size)

    flat_deriv = ni.gaussian_filter1d(flat_medianed, 1,
                                      order=1, axis=0)

    # min/max filter

    flat_max = ni.maximum_filter1d(flat_deriv, size=max_sep_order, axis=0)
    flat_min = ni.minimum_filter1d(flat_deriv, size=max_sep_order, axis=0)

    # mask for aperture boundray
    if pad is None:
        sl=slice()
    else:
        sl=slice(pad, -pad)

    flat_deriv_masked = np.zeros_like(flat_deriv)
    flat_deriv_masked[sl,sl] = flat_deriv[sl, sl]

    if flat_mask is not None:
        flat_deriv_pos_msk = (flat_deriv_masked > flat_max * 0.5) & flat_mask
        flat_deriv_neg_msk = (flat_deriv_masked < flat_min * 0.5) & flat_mask
    else:
        flat_deriv_pos_msk = (flat_deriv_masked > flat_max * 0.5)
        flat_deriv_neg_msk = (flat_deriv_masked < flat_min * 0.5)

    return dict(data=flat_deriv, #_bpix,
                pos_mask=flat_deriv_pos_msk,
                neg_mask=flat_deriv_neg_msk,
                )

def main():
    poly = load_file('monkey.ply')
    poly.BuildLinks(0)
    np = poly.GetNumberOfPoints()

    Ls = numpy.zeros((np, np))
    X, Y, Z = numpy_support.vtk_to_numpy(poly.GetPoints().GetData()).transpose()

    for i in xrange(np):
        for j in xrange(np):
            if i != j :
                if poly.IsEdge(i, j):
                    Ls[i, j] = -1
                    Ls[i, i] += 1

    Ls = numpy.matrix(Ls)
    Dx = X * Ls
    Dy = Y * Ls
    Dz = Z * Ls
    
    Dxg = gaussian_filter1d(Dx, 1)
    Dyg = gaussian_filter1d(Dy, 1)
    Dzg = gaussian_filter1d(Dz, 1)

    #Lsn = numpy.linalg.inv(Ls)

    #Xn = Lsn * Dxg
    #Yn = Lsn * Dyg
    #Zn = Lsn * Dzg

    X[:] = numpy.linalg.lstsq(Ls, Dxg.transpose())
    Y[:] = Dyg[0, :]
    Z[:] = Dzg[0, :]

    print X.shape
    
    V = numpy.zeros((np, 3))
    V[:,0] = X
    V[:,1] = Y
    V[:,2] = Z

    poly.GetPoints().SetData(numpy_support.numpy_to_vtk(V))
    poly.BuildLinks(0)
    poly.BuildCells()

    w = vtk.vtkSTLWriter()
    w.SetFileName('/tmp/teste.stl')
    w.SetInput(poly)
    w.Write()

def _erfimage_imshow_unified(bn, ch_idx, tmin, tmax, vmin, vmax, ylim=None,
                             data=None, epochs=None, sigma=None, order=None,
                             scalings=None, vline=None, x_label=None,
                             y_label=None, colorbar=False, cmap='RdBu_r',
                             vlim_array=None):
    """Plot erfimage topography using a single axis."""
    from scipy import ndimage
    _compute_ax_scalings(bn, (tmin, tmax), (0, len(epochs.events)))
    ax = bn.ax
    data_lines = bn.data_lines
    extent = (bn.x_t + bn.x_s * tmin, bn.x_t + bn.x_s * tmax, bn.y_t,
              bn.y_t + bn.y_s * len(epochs.events))
    this_data = data[:, ch_idx, :]
    vmin, vmax = (None, None) if vlim_array is None else vlim_array[ch_idx]

    if callable(order):
        order = order(epochs.times, this_data)

    if order is not None:
        this_data = this_data[order]

    if sigma > 0.:
        this_data = ndimage.gaussian_filter1d(this_data, sigma=sigma, axis=0)

    data_lines.append(ax.imshow(this_data, extent=extent, aspect='auto',
                                origin='lower', vmin=vmin, vmax=vmax,
                                picker=True, cmap=cmap,
                                interpolation='nearest'))