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

    mode1 = 'reflect'
    mode2 = ['reflect', 'reflect']

    assert_equal(sndi.gaussian_filter(arr, 1, mode=mode1),
                 sndi.gaussian_filter(arr, 1, mode=mode2))
    assert_equal(sndi.prewitt(arr, mode=mode1),
                 sndi.prewitt(arr, mode=mode2))
    assert_equal(sndi.sobel(arr, mode=mode1),
                 sndi.sobel(arr, mode=mode2))
    assert_equal(sndi.laplace(arr, mode=mode1),
                 sndi.laplace(arr, mode=mode2))
    assert_equal(sndi.gaussian_laplace(arr, 1, mode=mode1),
                 sndi.gaussian_laplace(arr, 1, mode=mode2))
    assert_equal(sndi.maximum_filter(arr, size=5, mode=mode1),
                 sndi.maximum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.minimum_filter(arr, size=5, mode=mode1),
                 sndi.minimum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.gaussian_gradient_magnitude(arr, 1, mode=mode1),
                 sndi.gaussian_gradient_magnitude(arr, 1, mode=mode2))
    assert_equal(sndi.uniform_filter(arr, 5, mode=mode1),
                 sndi.uniform_filter(arr, 5, mode=mode2))

    def gaussian_kernel(self,xvalues,yvalues,r200,normalization=100,scale=10,xres=200,yres=220,xmax=6.0,ymax=5000.0,adj=20):
        """
        Uses a 2D gaussian kernel to estimate the density of the phase space.
        As of now, the maximum radius extends to 6Mpc and the maximum velocity allowed is 5000km/s
        The "q" parameter is termed "scale" here which we have set to 10 as default, but can go as high as 50.
        "normalization" is simply H0
        "x/yres" can be any value, but are recommended to be above 150
        "adj" is a custom value and changes the size of uniform filters when used (not normally needed)
        """
        self.x_scale = xvalues/xmax*xres
        self.y_scale = ((yvalues+ymax)/(normalization*scale))/((ymax*2.0)/(normalization*scale))*yres

        img = np.zeros((xres+1,yres+1))
        self.x_range = np.linspace(0,xmax,xres+1)
        self.y_range = np.linspace(-ymax,ymax,yres+1) 

        for j in range(xvalues.size):
            img[self.x_scale[j],self.y_scale[j]] += 1
        
        #Estimate kernel sizes
        #Uniform
        #self.ksize = 3.12/(xvalues.size)**(1/6.0)*((np.var(self.x_scale[xvalues<r200])+np.var(self.y_scale[xvalues<r200]))/2.0)**0.5/adj
        #if self.ksize < 3.5:
        #    self.ksize = 3.5
        #Gaussian
        self.ksize_x = (4.0/(3.0*xvalues.size))**(1/5.0)*np.std(self.x_scale[xvalues<r200])
        self.ksize_y = (4.0/(3.0*yvalues.size))**(1/5.0)*np.std(self.y_scale[xvalues<r200])
        
        #smooth with estimated kernel sizes
        #img = ndi.uniform_filter(img, (self.ksize,self.ksize))#,mode='reflect')
        self.img = ndi.gaussian_filter(img, (self.ksize_y,self.ksize_x),mode='reflect')
        self.img_grad = ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x))
        self.img_inf = ndi.gaussian_gradient_magnitude(ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x)), (self.ksize_y,self.ksize_x))

def gborders(img, alpha=1.0, sigma=1.0):
    """Stopping criterion for image borders."""
    # The norm of the gradient.
    gausgradm = np.zeros(img.shape, np.double)
    gaussian_gradient_magnitude(img, sigma, output=gausgradm, mode='constant')
	
    thr = 2500
    mask = gausgradm > thr
    gausgradm[mask] = thr
    return 1.0/np.sqrt(1.0 + alpha*gausgradm)

    def pruning(self, skeleton_img, sigma):
        skeleton_img_mat = image_conversion.cv2array(skeleton_img)
        # Ausgabe-Array fuer das Ergebnis der Gradientberechnung
        gradient_output = numpy.empty_like(skeleton_img_mat)
        # Gradienten-Berechnung
        ndimage.gaussian_gradient_magnitude(skeleton_img_mat, sigma, gradient_output)
        # Normalisierung
        gradient_output /= gradient_output.max()
        # Array ins Bild umwandeln
        grad_img = image_conversion.array2cv(gradient_output)
        # Schwellwertbasierte Segmentierung des Gradientbildes
        dist_gradient_thresh = cv.CreateImage(cv.GetSize(grad_img), 8, 1)
        cv.InRangeS(grad_img, 0.6, 1, dist_gradient_thresh)

        return dist_gradient_thresh

def inverse_gaussian_gradient(image, alpha=100.0, sigma=5.0):
    """Inverse of gradient magnitude.

    Compute the magnitude of the gradients in the image and then inverts the
    result in the range [0, 1]. Flat areas are assigned values close to 1,
    while areas close to borders are assigned values close to 0.

    This function or a similar one defined by the user should be applied over
    the image as a preprocessing step before calling
    `morphological_geodesic_active_contour`.

    Parameters
    ----------
    image : (M, N) or (L, M, N) array
        Grayscale image or volume.
    alpha : float, optional
        Controls the steepness of the inversion. A larger value will make the
        transition between the flat areas and border areas steeper in the
        resulting array.
    sigma : float, optional
        Standard deviation of the Gaussian filter applied over the image.

    Returns
    -------
    gimage : (M, N) or (L, M, N) array
        Preprocessed image (or volume) suitable for
        `morphological_geodesic_active_contour`.
    """
    gradnorm = ndi.gaussian_gradient_magnitude(image, sigma, mode='nearest')
    return 1.0 / np.sqrt(1.0 + alpha * gradnorm)

def VesicleEdge_phc(img, x0, y0, r0, N=100, phi1=0, phi2=2 * np.pi, sigma=1):
    Xedge = np.empty(N)
    Yedge = np.empty(N)
    for i, phi in enumerate(np.linspace(phi1, phi2, N)):
        x = x0 + r0 * np.cos(phi)
        y = y0 + r0 * np.sin(phi)
        if x < 0:
            x = 0
            y = y0 + (x - x0) * np.tan(phi)
        elif x > img.shape[1] - 1:
            x = img.shape[1] - 1
            y = y0 + (x - x0) * np.tan(phi)
        if y < 0:
            y = 0
            x = x0 + (y - y0) / np.tan(phi)
        elif y > img.shape[0] - 1:
            y = img.shape[1] - 1
            x = x0 + (y - y0) / np.tan(phi)

        point1 = np.asarray(((y0, x0), (PIX_ERR, PIX_ERR)))
        point2 = np.asarray(((y, x), (PIX_ERR, PIX_ERR)))
        metric, metric_err, line = section_profile(img, point1, point2)
        grad = gaussian_gradient_magnitude(line, sigma)
        pos = np.argmax(grad)
        Xedge[i] = x0 + pos * np.cos(phi) * metric
        Yedge[i] = y0 + pos * np.sin(phi) * metric
    return Xedge, Yedge

 def remove_dust(self):
     self.seq = []
     for i in xrange(0,self.n_frames):
         # greyscale conversion
         img = np.average(self.reader.get_data(i),axis=2)
         self.seq.append(img)
     self.seq = np.array(self.seq)
     #var = np.var(self.seq, axis=0)
     #min = np.min(self.seq, axis=0)
     #max = np.max(self.seq, axis=0)
     #delta = max - min
     #var = stats.variation(self.seq, axis=0)
     #gmean = stats.gmean(self.seq, axis=0)
     a = np.average(self.seq, axis=0)
     #grad = ndimage.gaussian_gradient_magnitude(a , 0.25)
     #map = ndimage.prewitt(a)
     map = ndimage.gaussian_laplace(a,2.5) * ndimage.gaussian_gradient_magnitude(a , 0.25)
     cutoff = np.percentile(map,99.9)
     map[map<cutoff]=0
     map[map>0]=1
     #map = grad
     #map[map>300]=300
     fig = plt.figure(figsize=(20,8), frameon=False)
     fig.subplots_adjust(hspace=0)
     fig.subplots_adjust(wspace=0)
     ax1 = fig.add_subplot(1, 2, 1)
     ax1.imshow(map,interpolation='nearest')
     ax1.set_title('variance')
     ax2 = fig.add_subplot(1, 2, 2)
     ax2.imshow(self.seq[0], cmap='Greys_r',interpolation='nearest')
     ax2.set_title('img')
     fig.set_tight_layout(True)
     plt.show()

 def __init__(self, *args,**kwargs):
     im = kwargs.pop('image')
     super(Variance_no_D_Cmap,self).__init__(*args, **kwargs)
     
     fim = ndimage.gaussian_gradient_magnitude(im,2)
     fim = fim/np.std(fim)
     alpha = 1e0
     self.cmap = np.exp(-fim.flat[self.anchor['imask']]*alpha) + 1e-10

	def gaussian_kernel(self,xvalues,yvalues,r200,normalization=100,scale=10,res=200,adj=20,see=False):
		yres = 220
		#x_scale = (xvalues-np.min(xvalues))/np.max(xvalues-np.min(xvalues))*res
		#y_scale = ((yvalues-np.min(yvalues))/(normalization*scale))/np.max(xvalues-np.min(xvalues))*res
		self.x_scale = xvalues/6.0*res
		self.y_scale = ((yvalues+5000)/(normalization*scale))/(10000.0/(normalization*scale))*yres

		#img = np.zeros((int(np.max(x_scale))+1,int(np.max(y_scale))+1))
		img = np.zeros((res+1,yres+1))
		#x_range = np.linspace(np.min(xvalues),np.max(xvalues),int(np.max(x_scale))+1)
		#y_range = np.linspace(np.min(yvalues),np.max(yvalues),int(np.max(y_scale))+1)
		x_range = np.linspace(0,6,res+1)
		y_range = np.linspace(-5000,5000,yres+1) 

		for j in range(xvalues.size):
			img[self.x_scale[j],self.y_scale[j]] += 1

		#pcolormesh(img.T)
		#find ksize
		#xval = xvalues[np.where((xvalues<3) & (yvalues<2000) & (yvalues > -2000))]
		#yval = yvalues[np.where((xvalues<3) & (yvalues<2000) & (yvalues > -2000))]
		#x_scale2 = (xval-np.min(xval))/np.max(xval-np.min(xval))*res
		#y_scale2 = ((yval-np.min(yval))/(normalization*scale))/np.max(xval-np.min(xval))*res
		#xksize = 3.12/(xvalues.size)**(1.0/6.0)*((np.var(x_scale))/2.0)**0.5/adj
		#yksize = 3.12/(xvalues.size)**(1.0/6.0)*((np.var(y_scale))/2.0)**0.5/adj
	
		self.ksize = 3.12/(xvalues.size)**(1/6.0)*((np.var(self.x_scale[xvalues<r200])+np.var(self.y_scale[xvalues<r200]))/2.0)**0.5/adj
                self.ksize_x = (4.0/(3.0*xvalues.size))**(1/5.0)*np.std(self.x_scale[xvalues<r200])
                self.ksize_y = (4.0/(3.0*yvalues.size))**(1/5.0)*np.std(self.y_scale[xvalues<r200])
                if self.ksize < 3.5:
                    self.ksize = 3.5
		#ksize = 6.77588630223
		#print 'kernel size',ksize
		#img = ndi.uniform_filter(img, (self.ksize,self.ksize))#,mode='reflect')
                img = ndi.gaussian_filter(img, (self.ksize_y,self.ksize_x))#,mode='reflect')
                img_grad = ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x))
                img_inf = ndi.gaussian_gradient_magnitude(ndi.gaussian_gradient_magnitude(img, (self.ksize_y,self.ksize_x)), (self.ksize_y,self.ksize_x))
		
#		if see == True:
#			s = figure()
#			ax = s.add_subplot(111)
#			ax.pcolormesh(x_range,y_range,img.T)
#			show()
		return (x_range,y_range,img,np.abs(img_grad),np.abs(img_inf))

	def onEdgeFinding(self,evt):	
		if not self.panel.selectiontool.isTargeting('Auto Create Contours'):
			return 
		point = self.panel.view2image((evt.m_x,evt.m_y))

		ndarray = mrc.read(os.path.join(self.appionloop.params['rundir'], self.appionloop.imgtree[self.index]['filename']+'.dwn.mrc'))
		mrc.write(ndarray, os.path.join(self.appionloop.params['rundir'], 'beforefilter'+'.dwn.mrc'))

		negative = False	
		if self.filters:
			ndarray = ndimage.gaussian_filter(ndarray,1)
			ndarray = ndimage.gaussian_gradient_magnitude(ndarray,2)
			markers = []
			for i in range(3):
				for j in range(3):
					if i!=0 or j!=0:
						markers.append((point[0]-1+i,point[1]-1+j))
			markers = (1,2,3,4,5,6,7,8)
			#ndarray = ndimage.watershed_ift(ndarray,markers)
			ndarray = ndimage.laplace(ndarray)
			ndarray = ndimage.gaussian_filter(ndarray,1)
			#ndarray = apImage.preProcessImage(ndarray,params=self.appionloop.params)
			negative = True 
		mrc.write(ndarray, os.path.join(self.appionloop.params['rundir'], 'afterfilter'+'.dwn.mrc'))

		delta = .1
		targets = []
		radius = 20
		size = 50
		rangeSize = 50
		maker = PixelCurveMaker()
		maker._init_(size,rangeSize);
		for theta in range(size):
			theta +=0
			theta*=math.pi*2/rangeSize
			for rad in range(size):
				try:
					if negative:
						maker.addData(theta,rad,127-ndarray[int(point[1]+rad*math.sin(theta))][int(point[0]+rad*math.cos(theta))])
					else:
						maker.addData(theta,rad,ndarray[int(point[1]+rad*math.sin(theta))][int(point[0]+rad*math.cos(theta))])
				except IndexError:
					maker.addData(theta,rad,0)
		maker.makeCalculations()
		s = self.filterSelectorChoices[self.filterSelector.GetSelection()]	
		dilate = 2
		if s == 'Latex Bead':
			dilate = 0
		for theta in range(size):
			theta += 0
			theta*=math.pi*2/rangeSize
			targets.append((point[0]+(dilate+maker.getData(theta))*math.cos(theta),point[1]+(dilate+maker.getData(theta))*math.sin(theta)))
		self.addPolyParticle(targets)
		#this section draws all of the contours that the algorithm considers - useful for debugging
		'''

def getGradientVideo(I, IDims, sigma = 1):
    GV = np.zeros(I.shape)
    for i in range(I.shape[0]):
        X = np.reshape(I[i, :], IDims)
        G = rgb2gray(X, False)
        GM = gaussian_gradient_magnitude(G, sigma)
        F = np.zeros(IDims)
        for k in range(F.shape[2]):
            F[:, :, k] = GM
        GV[i, :] = F.flatten()
    return GV

def hipass(image,gaussrad,dilrat):
    # 
    # A derivative of gaussian (suggested radii of 0.5) to get edges
    # followed by a morphological dilate to widen the effect and pick up extra pixels on edges
    # This can be used as a mask for 

    timg = np.copy(image).reshape(96,96)
    
    gradim = ndimage.gaussian_gradient_magnitude(timg,gaussrad)
    
    digradim= ndimage.morphology.grey_dilation(gradim,size=(dilrat,dilrat))
    
    return digradim.flatten()

def test_multiple_modes_gaussian_gradient_magnitude():
    # Test gaussian_gradient_magnitude filter for multiple
    # extrapolation modes
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    expected = np.array([[0.04928965, 0.09745625, 0.06405368],
                         [0.23056905, 0.14025305, 0.04550846],
                         [0.19894369, 0.14950060, 0.06796850]])

    modes = ['reflect', 'wrap']

    calculated = sndi.gaussian_gradient_magnitude(arr, 1, mode=modes)

    assert_almost_equal(expected, calculated)

def task2_4():
    img = normalize_intensity(imread(CAMERAMAN))
    img = img[30:95, [i for i in range(80, 160)]]  # select subsection of image
    vel_x, vel_y = gradient(f=img)

    magn_img = gaussian_gradient_magnitude(img, 3)
    output_path = os.path.join(OUTPUT_DIR, "2_4_gradien_magnitude_" + os.path.split(CAMERAMAN)[-1])
    imsave(output_path, magn_img)

    dim_x, dim_y = len(img[0]), len(img)
    x, y = range(dim_x), range(dim_y)
    x, y = meshgrid(x, y)
    plt.figure()
    imgplot = plt.imshow(img)
    imgplot.set_cmap('gray')
    plt.ylim(dim_y, 0)
    plt.quiver(x, y, vel_x, vel_y, pivot='middle')
    plt.show()

def wall_points_pix_3(img, refsx, axisy, sigma):
    N = 2
    refs = np.array([])
    axisy = np.array([113, 115])
    for i, refx in enumerate(refsx):
        prof = img[:, refx]
        mid = axisy[i]
        #        mid = len(prof)/2
        filtered = ndimage.gaussian_gradient_magnitude(ndimage.sobel(prof), sigma)

        refy = np.asarray((np.argmax(filtered[:mid]), np.argmax(filtered[mid:]) + mid))
        dref = np.asarray([PIX_ERR, 0])
        rx = np.tile(refx, N)
        xy = np.column_stack((refy, rx))  # .flatten()
        drefe = np.repeat(np.expand_dims(dref, 0), N, 0)
        ref = np.concatenate((xy, drefe), 1)
        refs = np.append(refs, ref)
    return refs.reshape(-1, 2, 2)

def fitFirstCTFNode(pow, rpixelsize, defocus, ht):
	filter = ndimage.gaussian_filter(pow,3)
	grad = ndimage.gaussian_gradient_magnitude(filter,3)
	thr = imagefun.threshold(grad,grad.mean()+3*grad.std())
	if defocus:
		z = abs(defocus)
		s = calculateFirstNode(ht,z)
		dmean = max(0.8*s/rpixelsize, 30)
	else:
		shape = pow.shape
		r = 20
		center = ( shape[0] / 2, shape[1] / 2 ) 
		grad[center[0]-r: center[0]+r, center[1]-r: center[1]+r] = 0
		peak = ndimage.maximum_position(grad)
		dmean = math.hypot(peak[0] - center[0], peak[1] - center[1])
	drange = max(dmean / 4, 10)
	eparams = find_ast_ellipse(grad,thr,dmean,drange)
	if eparams:
		z0, zast, ast_ratio, alpha = getAstigmaticDefocii(eparams,rpixelsize, ht)
		return z0,zast,ast_ratio, alpha, eparams

	def predict_proba_image(self, array_data, w_x, w_y):
		"""Predict class probabilities for X"""
		all_proba = []

		if(self.use_geodesic):
			pan = array_data[0]
			self.geodesic_cost = nd.gaussian_gradient_magnitude(pan, self.geodesic_sigma)


		if(self.n_steps_simple is None and self.n_steps_proba is None):
			for i in range(self.n_forests):
				if ((i != 0) and self.add_previous_prob):
					if(self.use_geodesic):
						proba = self.geodesic(proba)
					array_data = numpy.concatenate((array_data,proba))
				proba = self.forests_[i].predict_proba_image(array_data, w_x, w_y)
				all_proba.append(proba)
		else:
			i = 0
			done = True
			for step_proba in range(self.n_steps_proba):
				for step_simple in range(self.n_steps_simple):
					if (step_proba != 0) and (step_simple == 0):
						proba = self.forests_[i].predict_proba_image(array_data, w_x, w_y)
						if(self.use_geodesic):
							proba = self.geodesic(proba)
						#if use_geodesic
						array_data = numpy.concatenate((array_data,proba))
					#if (step_proba != 0) and (step_simple=0):
					proba = self.forests_[i].predict_proba_image(array_data, w_x, w_y)
					i +=1
				#for step_simple
			#for step_proba

		if(self.fusion == "mean"):
			for j in range(1, len(all_proba)):
				proba += all_proba[j]
			return proba / self.n_forests
		else: # (fusion =="last"):
			return proba

def gdt(binaryImage, image, voxelSpacing=(1.0, 1.0, 1.0), subtractInside=True, includeEDT=True,
        gamma=32, numIterations=3, boundingMask=None, magnitudes=None):
    """
        Runs geodesic distance transform based on gradient magnitues as the geodesic cost
        @param boundingMask is a mask provided where its bounding box will be used to restrict
               where the distance transform is run for performance gain
    """
    if numpy.count_nonzero(binaryImage) == 0:
        print "[WARN]--->No binary data provided to GDT"
        distanceMap = numpy.zeros(binaryImage.shape)
        distanceMap[:] = numpy.NaN
        return distanceMap

    if magnitudes is None:
        image = numpy.asarray(image, dtype=numpy.float64)
        magnitudes = gaussian_gradient_magnitude(image, 0.5)
        if includeEDT and gamma != 1:
            magnitudes *= gamma
    minBounds, maxBounds = None, None
    if boundingMask is not None:
        minBounds, maxBounds = get_bounding_box(boundingMask)
        binaryImage = binaryImage[minBounds[0]:maxBounds[0], minBounds[1]:maxBounds[1], minBounds[2]:maxBounds[2]]
        magnitudes = magnitudes[minBounds[0]:maxBounds[0], minBounds[1]:maxBounds[1], minBounds[2]:maxBounds[2]]

    distanceMap = geodesic_distance_transform(binaryImage, magnitudes,
                                              numIterations=numIterations, spacing=voxelSpacing, includeEDT=includeEDT)
    if subtractInside:
        distanceMap -= geodesic_distance_transform(logical_not(binaryImage), magnitudes,
                                                   numIterations=numIterations, spacing=voxelSpacing,
                                                   includeEDT=includeEDT)
    if magnitudes is None and includeEDT and gamma != 1:
        # make it more comparable when different gamma are used
        distanceMap /= gamma

    if boundingMask is not None:
        temp = numpy.zeros(boundingMask.shape)
        temp[minBounds[0]:maxBounds[0], minBounds[1]:maxBounds[1], minBounds[2]:maxBounds[2]] = distanceMap
        distanceMap = temp

    return distanceMap

def gdt( img, mask, includeEDT=True, l=1.0 ):
    if mask.sum() == 0:
        return irtk.zeros(img.get_header())
    voxelSpacing = img.header['pixelSize'][:3][::-1]
    grad =  irtk.Image( nd.gaussian_gradient_magnitude(img, 0.5),
                        img.get_header() )
    #irtk.imwrite("gradBefore.nii.gz",grad)
    grad = l*grad.saturate().rescale(0.0,1.0).as3D()
    #irtk.imwrite("gradAfter.nii.gz",grad)
    # distanceMap = geodesic_distance_transform( mask,
    #                                            grad,
    #                                            numIterations=3,
    #                                            spacing=voxelSpacing,
    #                                            includeEDT=includeEDT )
    # distanceMap -= geodesic_distance_transform( logical_not(mask),
    #                                             grad,
    #                                             numIterations=3,
    #                                             spacing=voxelSpacing,
    #                                             includeEDT=includeEDT )
    # return irtk.Image(distanceMap,img.get_header())

    # distanceMaps = Parallel(n_jobs=-1)(delayed(_geodesic_distance_transform)( m,
    #                                                                          grad,
    #                                                                          numIterations=3,
    #                                                                          spacing=voxelSpacing,
    #                                                                          includeEDT=includeEDT )
    #                                                              for m in [mask,
    #                                                                        logical_not(mask)]
    #                                                                          )
    # res = irtk.Image(distanceMaps[0]-distanceMaps[1],img.get_header())

    res = irtk.Image( _geodesic_distance_transform( mask,
                                                    grad,
                                                    numIterations=3,
                                                    spacing=voxelSpacing,
                                                    includeEDT=includeEDT ),
                      img.get_header() ).as3D()

    return res

def wall_points_pix_2(img, refsx, sigma):
    """
    
    @param img:
    @param refsx:
    """
    N = 2  # number of walls to find
    refs = np.array([])
    for refx in refsx:
        prof = img[:, refx]
        gradprof = ndimage.gaussian_gradient_magnitude(prof, sigma)
        start, end = split_two_peaks(gradprof, 1)
        if start > end:
            start, end = end, start
        refy = split_two_peaks(prof[start:end], -1) + start
        dref = np.asarray([PIX_ERR, 0])
        rx = np.tile(refx, N)
        xy = np.column_stack((refy, rx))  # .flatten()
        drefe = np.repeat(np.expand_dims(dref, 0), N, 0)
        ref = np.concatenate((xy, drefe), 1)
        refs = np.append(refs, ref)
    return refs.reshape(-1, 2, 2)