def use_ridge_regression():
    # Read the command line arguments.
    args = process_command_line()
    im_orig = numpy.squeeze(vigra.readImage("cc_90.png"))
    im = kernel_ridge_regression(im_orig, args.tau, args.rho, args.gamma)
    vigra.impex.writeImage(im, "res.png")
    im_true = numpy.squeeze(vigra.readImage("charlie-chaplin.jpg"))
    print "MSE: ", calc_mse(im, im_true)

    return 0

def importStack(path,fname):
	absname = path +fname
	zsize = vigra.impex.numberImages(absname)
	im =vigra.readImage(absname, index = 0, dtype='FLOAT')
	#vol = np.zeros([im.height,im.width,zsize])
	vol = np.memmap('tmpVolDat/' + fname[0:-4],dtype='float64',mode = 'w+', shape = (im.height,im.width,zsize))
	#raise('hallo')
	for i in range(zsize):
		print("importing slice " + str(i) + ' of file '+fname)
		im=np.squeeze(vigra.readImage(absname, index = i, dtype='FLOAT'))
		vol[:,:,i] = im
	vol.flush()
	return vol

def bayesian_hp_optimization(Q):
    im = numpy.squeeze(vigra.readImage("cc_90.png"))
    im_orig = numpy.squeeze(vigra.readImage("charlie-chaplin.jpg"))
    P = []
    E = []
    # initialize
    f = open("cache/bayes_opt.txt","w")
    start = time.time()
    #for i in range(20):
    #    interpolation = kernel_ridge_regression( im, Q[i,2], Q[i,0], Q[i,1] )
    #    P.append( [Q[i,0], Q[i,1], Q[i,2]] )
    #    E.append( calc_mse(im_orig, interpolation) )
    #    # save result
    #    res = str(Q[i,0]) + str(" ") + str(Q[i,1]) + str(" ") + str(Q[i,2]) + str(" ") + str(E[i]) + '\n'
    #    f.write(res)
    #    f.flush()

    save = numpy.loadtxt("cache/save.txt")
    for d in save:
        P.append( [d[0], d[1], d[2]] )
        E.append( d[3] )

    # TODO should we remove known vals from Q ?
    # remove known values from Q
    # Q = numpy.delete(Q, numpy.arange(20), axis=0)
    # parameter for the matern regression
    sig_rho     = 8. / 10
    sig_gamma   = 3. / 10
    sig_tau     = .9 / 10
    lambd       = .3
    for i in range(20):
        mse, var = matern_regression(Q, P, E, sig_rho, sig_gamma, sig_tau, lambd)
        #utility = numpy.divide( mse, numpy.sqrt(var) )
        utility = numpy.abs(numpy.divide( mse, var ) )
        best_hp = numpy.nanargmin(utility)
        P.append( Q[best_hp])
        print Q[best_hp]
        print utility[best_hp], mse[best_hp], var[best_hp]
        interpolation = kernel_ridge_regression( im, Q[best_hp,2], Q[best_hp,0], Q[best_hp,1])
        E.append( calc_mse(im_orig, interpolation))
        res = str(Q[best_hp,0]) + str(" ") + str(Q[best_hp,1]) + str(" ") + str(Q[best_hp,2]) + str(" ") + str(E[-1]) + '\n'
        f.write(res)
        f.flush()
    stop = time.time()
    print "Bayesian parameter optimization took %.02f seconds." % (stop-start)
    best_hp = numpy.argmin(E)
    f.close()
    return P[best_hp], E[best_hp]

    def load_image(self, filename=None):
        # if filename is None:
        #    filename = self.image_name

        # get input image
        self.image_name = filename
        self.set_image(vigra.readImage(filename))

 def setImage(self, imagefilename, retainView=False):
     self.image = vigra.readImage(imagefilename)
     shapefactor = self.image.shape[0]/5000
     self.imagedisplay = vigra.sampling.resizeImageNoInterpolation(self.image, (self.image.shape[0]/shapefactor, self.image.shape[1]/shapefactor))
     self.imagedisplay = vigra.colors.brightness(vigra.colors.linearRangeMapping(self.imagedisplay), 35.)
     super(S57QImageViewer, self).setImage(self.imagedisplay.qimage(), retainView)
     self.geoimage = GeoImage(imagefilename)

def kernel_ridge_regression(tau, sigma):
    # Load the image.
    im_orig = numpy.squeeze(vigra.readImage("cc_90.png"))

    # Make a copy, so both the original and the regressed image can be shown afterwards.
    im = numpy.array(im_orig)

    # Find the known pixels and the pixels that shall be predicted.
    known_ind = numpy.where(im != 0)
    unknown_ind = numpy.where(im >= 0)
    known_x = numpy.array(known_ind).transpose()
    known_y = numpy.array(im[known_ind])
    pred_x = numpy.array(unknown_ind).transpose()

    # Train and predict with the given regressor.
    start = time.time()
    print "training..."
    r = KernelRidgeRegressor(tau, sigma)
    r.train(known_x, known_y)
    print "done training"
    # pickle.dump(r, open("regressor.p", "wb"))
    # r = pickle.load(open("regressor.p", "rb"))
    print "predicting..."
    pred_y = r.predict(pred_x)
    print "done predicting"

    # Write the predicted values back into the image and show the result.
    im[unknown_ind] = pred_y
    stop = time.time()
    print "Train and predict took %.02f seconds." % (stop-start)
    vigra.impex.writeImage(im, "res.png")

 def load_image(self, file_):
     if not self._image_cache.has_key(file_):
         image = vigra.readImage(file_)
         # numpy array convention
         image.swapaxes(0, 1)
         self._image_cache[file_] = np.squeeze(image)
     return self._image_cache[file_]

    def addBackgroundWithFrame(self, bgImageFilename, container = True, **params):
        """fe.addBackgroundWithFrame(bgImageFilename, depth = 85, ...)

        Adds a picture object to the fig.File, framed by an additional
        rectangle.  See addROIRect() and addImage().  If no roi is
        given (via a keyword parameter), the image file is opened
        (using readImage) and its size is used to initialize a
        BoundingBox positioned at the origin.

        Returns the pair (bgImage, bgRect) of both added fig objects."""

        if not params.has_key("roi") and not self.roi:
            size = readImage(bgImageFilename).size()
            params["roi"] = Rect2D(size)
        if container == True:
            container = self.f

        if not params.has_key("depth"):
            params["depth"] = 1

        bgRect = self.addROIRect(**params)

        bgImage = fig.PictureBBox(0, 0, 1, 1, bgImageFilename)
        bgImage.points = list(bgRect.points)
        bgImage.depth = 999
        container.append(bgImage)

        return bgImage, bgRect

def extractMarkedNodes(filename):
    # extracts markers from the raw images
    # returns a dict, with color as key and a list
    # of marker center coords as value

    print "processing file", filename
    im = vigra.readImage(filename)
    colored1 = im[..., 0] != im[..., 1]
    colored2 = im[..., 1] != im[..., 2]
    colored3 = im[..., 2] != im[..., 0]
    colored = numpy.logical_or(colored1, colored2)
    colored = numpy.logical_or(colored, colored3)
    cc = vigra.analysis.labelImageWithBackground(colored.astype(numpy.uint8))
    # take the center pixel for each colored square
    feats = vigra.analysis.extractRegionFeatures(colored.astype(numpy.float32), cc, ["RegionCenter"])
    center_coords = feats["RegionCenter"][1:][:].astype(numpy.uint32)
    center_coords_list = [center_coords[:, 0], center_coords[:, 1]]
    im_centers = numpy.asarray(im[center_coords_list])
    # print im_centers
    # struct = im_centers.view(dtype='f4, f4, f4')

    # colors, indices = numpy.unique(struct, return_inverse=True)
    # print colors, indices, colors.shape
    centers_by_color = {}
    for iindex in range(center_coords.shape[0]):
        center = (center_coords[iindex][0], center_coords[iindex][1])
        #print center, index
        color = colors[tuple(im_centers[iindex].astype(numpy.uint8))]
        #centers_by_color.setdefault(tuple(im_centers[iindex]), []).append(center)
        centers_by_color.setdefault(color, []).append(center)
        
    print centers_by_color
    return centers_by_color

    def get_data(self, data_nr):
        """Returns the dataset.

        :param data_nr: number of dataset
        :return: the dataset
        """
        if self._datatype(data_nr) == "hdf5" or self.is_internal(data_nr):
            return vigra.readHDF5(self.get_data_path(data_nr), self.get_data_key(data_nr))
        else:
            return vigra.readImage(self.get_data_path(data_nr))

 def _execute_5d(self, roi, result):
     t_start, x_start, y_start, z_start, c_start = roi.start
     t_stop, x_stop, y_stop, z_stop, c_stop = roi.stop
     
     for result_t, t in enumerate( range( t_start, t_stop ) ):
         file_name = self.fileNameList[t]
         for result_z, z in enumerate( range( z_start, z_stop ) ):
             img = vigra.readImage( file_name, index=z )
             result[result_t, :, :, result_z, :] = img[ x_start:x_stop,
                                                        y_start:y_stop,
                                                        c_start:c_stop ]
     return result

def random_hp_optimization(Q):
    f = open("cache/rand_opt.txt","w")
    image = numpy.squeeze(vigra.readImage("cc_90.png"))
    im_orig = numpy.squeeze(vigra.readImage("charlie-chaplin.jpg"))
    start = time.time()
    P = []
    E = []
    N = Q.shape[0]
    rand_indices = numpy.random.randint(0, N, size = 40)
    for i in rand_indices:
        interpolation = kernel_ridge_regression( image, Q[i,2], Q[i,0], Q[i,1] )
        P.append( [Q[i,0], Q[i,1], Q[i,2]] )
        E.append( calc_mse(im_orig, interpolation) )
        res = str(Q[i,0]) + str(" ") + str(Q[i,1]) + str(" ") + str(Q[i,2]) + str(" ") + str(E[-1]) + '\n'
        f.write(res)
        f.flush()
    rand_hp = numpy.argmin(E)
    stop = time.time()
    print "Random parameter optimization took %.02f seconds." % (stop-start)
    f.close()
    return P[rand_hp], E[rand_hp]

def main(filename, biScale = 1.6, saddleThreshold = 0.2):
    """Creates an initial GeoMap ("level 0" of irregular pyramid) using
    subpixel watersheds on a Gaussian gradient boundary indicator and
    creates a Workspace from that."""

    import bi_utils
    img = vigra.readImage(filename)
    gm, grad = bi_utils.gaussianGradient(img, biScale)
    wsm = maputils.subpixelWatershedMap(
        gm, saddleThreshold = saddleThreshold)

    return Workspace(wsm, img, bi = gm)

def volume_from_dir(dirpattern, output_filepath, offset=0, nfiles=None):
    filelist = glob.glob(dirpattern)
    filelist = sorted(filelist, key=str.lower) #mwahaha, 10000<9000
    begin = offset
    if nfiles is not None and offset+nfiles<len(filelist):
        end=offset+nfiles
    else:
        end = len(filelist)
    filelist = filelist[begin:end]
    nx, ny = vigra.readImage(filelist[0]).squeeze().shape
    dt = vigra.readImage(filelist[0]).dtype
    nz = len(filelist)
    volume = numpy.zeros((nx, ny, nz, 1), dtype=dt)
    
    for i in range(len(filelist)):
        volume[:, :, i, 0] = vigra.readImage(filelist[i]).squeeze()[:]
        
    outfile = h5py.File(output_filepath, "w")
    outfile.create_dataset("data", data=volume)
    outfile.close()
    return volume

def makeTif(pathSearch, pathSave, filename):  
    frameNum = os.path.splitext(filename)[0][-5:]; begName = os.path.splitext(filename)[0][:-5]; ext = os.path.splitext(filename)[1]
    newFrameNum = '{:0>5}'.format(int(frameNum)*30)
    try:
        im=vigra.readImage(os.path.join(pathSearch, filename))
    except IOError:
        print "File pbl with file ", filename
        return 0
    else:
        im2=vigra.Image(im, dtype=np.uint8)
        im2[im2>0]=1
        im2.writeImage(os.path.join(pathSave, 'Mask_'+begName+newFrameNum+ext))
        return 1

    def _execute_5d(self, roi, result):
        # roi is in tzyxc order.
        t_start, z_start, y_start, x_start, c_start = roi.start
        t_stop, z_stop, y_stop, x_stop, c_stop = roi.stop

        # Use *enumerated* range to get global t coords and result t coords
        for result_t, t in enumerate( range( t_start, t_stop ) ):
            file_name = self.fileNameList[t]
            for result_z, z in enumerate( range( z_start, z_stop ) ):
                img = vigra.readImage( file_name, index=z )
                result[result_t, result_z, :, :, :] = img[ x_start:x_stop,
                                                           y_start:y_stop,
                                                           c_start:c_stop ].withAxes( *'yxc' )
        return result

def dense_reconstruction_slice(probs_slice_path, skeletons_volume, rf):
    probs = vigra.readImage(probs_slice_path)
    probs = np.squeeze(probs)
    probs = np.array(probs) # get rid of axistags...
    # need to invert for watersheds
    probs = 1. - probs

    # Threshold and sigma hardcoded to values suited for the google pmaps
    seg_wsdt = watersheds_dt(probs, 0.15, 2.6)

    # this may happen for the black slices...
    if np.unique(seg_wsdt).shape[0] == 1:
        return np.zeros_like(probs, dtype = np.uint32)

    rag = vigra.graphs.regionAdjacencyGraph(
            vigra.graphs.gridGraph(seg_wsdt.shape[0:2]), seg_wsdt )

    # +1 because we dont have a zero in overseg
    merge_nodes = np.zeros(rag.nodeNum+1, dtype = np.uint32)

    gridGraphEdgeIndicator = vigra.graphs.implicitMeanEdgeMap(rag.baseGraph,
            probs)
    edge_feats = rag.accumulateEdgeStatistics(gridGraphEdgeIndicator)
    edge_feats = np.nan_to_num(edge_feats)
    edge_probs = rf.predict_proba(edge_feats)[:,1]
    probs_thresh = 0.5

    for skel_id in np.unique(skeletons_volume):
        if skel_id == 0:
            continue
        skel_c = np.where(skeletons_volume == skel_id)
        for i in xrange(skel_c[0].size):
            seg_id = seg_wsdt[ skel_c[0][i], skel_c[1][i] ]
            merge_nodes[seg_id] = skel_id
            root = rag.nodeFromId( long(seg_id) )
            nodes_for_merge = [root]
            already_merged  = [root]
            while nodes_for_merge:
                u = nodes_for_merge.pop()
                for v in rag.neighbourNodeIter(u):
                    edge = rag.findEdge(u,v)
                    #print edge_mean_probs[edge.id]
                    if edge_probs[edge.id] <= probs_thresh and not v.id in already_merged:
                        merge_nodes[v.id] = skel_id
                        nodes_for_merge.append(v)
                        already_merged.append(v.id)
    return rag.projectLabelsToBaseGraph(merge_nodes)

 def __call__(self, filename, out_path):
     if not os.path.exists(out_path):
         os.makedirs(out_path)
         print 'made %s' % out_path
         
     colorin = vigra.readImage(filename)
     filename_base, extension = os.path.splitext(os.path.basename(filename))
     
     col_dec = self.colorDeconv(colorin)
     channels = ['h', 'e', 'dab']
     for i in range(3):
         new_filename = os.path.join(out_path, 
                                     filename_base + '__%s' % channels[i] + extension)
         vigra.impex.writeImage(col_dec[:,:,i], new_filename) 
         print 'written %s' % new_filename
         
     return

def getImageByName(imagefilename, topleft=None, bottomright=None, linearrangemapping=True):
	"""
		Returns a defined crop of the file as vigra.Image. 
	"""
	base, ext = os.path.splitext(imagefilename)
	if not ext[1:] in getVigraExts():
		raise IOError('unsupported file extension %s' % ext[1:])
	image = vigra.readImage(imagefilename)
	if topleft == None and bottomright == None:
		return image
	# determine crop size
	x_size = bottomright[0] - topleft[0]
	y_size = bottomright[1] - topleft[1]
	# crop
	out = image[topleft[0]:topleft[0]+x_size, topleft[1]:topleft[1]+y_size].copy()
	if linearrangemapping:
		out = vigra.colors.linearRangeMapping(out)
	return out

def get_extensions(in_folder):

    max_width = 0
    max_height = 0    
    image_names = os.listdir(in_folder)
    image_names = sorted(filter(lambda x: os.path.splitext(x)[-1].lower() in ['.tif', '.tiff', '.png', '.jpg'], image_names))
    
    for image_name in image_names:
        img = vigra.readImage(os.path.join(in_folder, image_name))
        width = img.shape[0]
        height = img.shape[1]        
        print '%s: %i, %i' % (image_name, width, height)
        
        max_width = max(width, max_width)
        max_height = max(height, max_height)
        
    print 'maximal extensions: %i, %i' % (max_width, max_height)
    return max_width, max_height