def test_set_color_with_alpha():
    img = np.zeros((10, 10))

    rr, cc, alpha = line_aa(0, 0, 0, 30)
    set_color(img, (rr, cc), 1, alpha=alpha)

    # Wrong dimensionality color
    assert_raises(ValueError, set_color, img, (rr, cc), (255, 0, 0), alpha=alpha)

    img = np.zeros((10, 10, 3))

    rr, cc, alpha = line_aa(0, 0, 0, 30)
    set_color(img, (rr, cc), (1, 0, 0), alpha=alpha)

def test_line_equal_aliasing_horizontally_vertically():
    img0 = np.zeros((25, 25))
    img1 = np.zeros((25, 25))

    # Near-horizontal line
    rr, cc, val = line_aa(10, 2, 12, 20)
    img0[rr, cc] = val

    # Near-vertical (transpose of prior)
    rr, cc, val = line_aa(2, 10, 20, 12)
    img1[rr, cc] = val

    # Difference - should be zero
    assert_array_equal(img0, img1.T)

def test_line_aa_vertical():
    img = np.zeros((10, 10))

    rr, cc, val = line_aa(0, 0, 9, 0)
    img[rr, cc] = val

    img_ = np.zeros((10, 10))
    img_[:, 0] = 1

    assert_array_equal(img, img_)

def test_line_aa_horizontal():
    img = np.zeros((10, 10))

    rr, cc, val = line_aa(0, 0, 0, 9)
    img[rr, cc] = val

    img_ = np.zeros((10, 10))
    img_[0, :] = 1

    assert_array_equal(img, img_)

def test_line_aa_horizontal():
    img = np.zeros((10, 10))

    rr, cc, val = line_aa(0, 0, 0, 9)
    set_color(img, (rr, cc), 1, alpha=val)

    img_ = np.zeros((10, 10))
    img_[0, :] = 1

    assert_array_equal(img, img_)

 def crop_to_sequence(self, track_list, currSharp):
     '''
     In the crop filename, I add the number of the image in the galerie so that they're by default ordered in time
     
     I need the bounding box of the membrane and also to use whitefield images.
     '''
     
     for track in track_list:
         id_=track.id
         lstFrames=sorted(track.lstPoints.keys(), key=itemgetter(0))
         rr=[]; cc=[]; val=[]; nextCoord=None
         for i, el in enumerate(lstFrames):
             im, cell_id=el
             coordonnees=track.lstPoints[(im, cell_id)] if nextCoord==None else nextCoord
             try:
                 nextCoord=track.lstPoints[lstFrames[i+1]]
             except IndexError:
                 continue
             else:
                 r,c,v=draw.line_aa(coordonnees[0], coordonnees[1],nextCoord[0], nextCoord[1])
                 rr.extend(r); cc.extend(c); val.extend(v)
                 
         for im, cell_id in lstFrames:
             #renumbering according to xb screen/PCNA image numbering
             local_im=im+1
             
             #draw a dot on the cell which is followed
             cell_x, cell_y=track.lstPoints[(im, cell_id)]
             dot_rr, dot_cc=draw.circle(cell_x, cell_y, radius=3)
                 
             image_name= self.settings.imageFilename.format(self.well, local_im)
             image=vi.readImage(os.path.join(self.settings.allDataFolder, self.plate, 'analyzed', self.well, 'images/tertiary_contours_expanded', image_name))
             
             #X,x,x__, Y,y,y__=self._newImageSize(crop_coordinates)
             x__=self.settings.XMAX; y__=self.settings.YMAX; x=0; y=0; X=self.settings.XMAX; Y=self.settings.YMAX
             
             croppedImage = VigraArray((x__, y__, 3), dtype=np.dtype('float32'))
             croppedImage=image[x:X, y:Y]  
             croppedImage[rr,cc,0]=np.array(val)*255
             
     #If there is a sharp movement, the cell center is pinky red
             if im in currSharp[id_]:
                 croppedImage[dot_rr, dot_cc, 0]=242
                 croppedImage[dot_rr, dot_cc, 1]=21
                 croppedImage[dot_rr, dot_cc, 2]=58
             else:
     #If not, it is green
                 croppedImage[dot_rr, dot_cc, 1]=255
             
             vi.writeImage(croppedImage, \
                           os.path.join(self.outputFolder, self.plate, 'galerie',
                                        self.settings.outputImage.format(self.plate, self.well.split('_')[0],id_, im)),\
                           dtype=np.dtype('uint8'))
             
     return

def test_line_aa_diagonal():
    img = np.zeros((10, 10))

    rr, cc, val = line_aa(0, 0, 9, 6)
    img[rr, cc] = 1

    # Check that each pixel belonging to line,
    # also belongs to line_aa
    r, c = line(0, 0, 9, 6)
    for x, y in zip(r, c):
        assert_equal(img[r, c], 1)

def projection_on_magnet(ends_of_strawtubes, x_resolution=1., y_resolution=1., z_magnet=3070.):
    
    tubes = []
    
    for i in range(len(ends_of_strawtubes)):
        
        tubes.append(ends2params(ends_of_strawtubes[i]))

    lines = []
    z3 = z_magnet
    for i in range(len(tubes)):
        
        for j in range(i+1, len(tubes)):
        
            t1 = tubes[i]
            t2 = tubes[j]
            
            if (t1[2] != t2[2]) and (np.sum(t1[1] - t2[1]) < 0.01):
                
                start1 = t1[0]
                start2 = t2[0]
                z1 = t1[2]
                z2 = t2[2]
                start3 = (start2 - start1) / (z2 - z1) * (z3 - z2) + start2
                
                if (start3[1] > -500) and (start3[1] < 500) and (start3[0] < -225) and (start3[0] > -275):
                
                    lines.append([start3, t1[1]])
    
    x_len = int(600 * x_resolution)
    y_len = int(1200 * y_resolution)
    line_len = 500

    matrix = np.zeros([x_len, y_len])
    
    for line in lines:
        
        rr, cc, val = line_aa(int(round(line[0][0] * x_resolution)) + x_len / 2, int(round(line[0][1] * y_resolution)) + y_len / 2,\
                              int(round((line[0][0] + line[1][0] * line_len) * x_resolution)) + x_len / 2,\
                              int(round((line[0][1] + line[1][1] * line_len) * y_resolution)) + y_len / 2)
        matrix[rr, cc] += 1
        
    return matrix

def filter_contours_lines(v0, v1, lines, img, dist): #linie zaczynajace sie w danym punkcie w otoczeniu dist pikseli
  print "v0, v1:" , v0, v1
  max_line_length = 0
  for tmp_line in lines:
    tmp0, tmp1 = tmp_line
    tp0 = (tmp0[1], tmp0[0])
    tp1 = (tmp1[1], tmp1[0])
    if distance(tmp0, tmp1) > max_line_length:
      max_line_length = distance(tmp0, tmp1)
      print "NEW_MAX:", tmp0, tmp1
    print "tp0 tp1: " ,tp0, tp1
    print "v0 tp0", distance(v0, tp0)
    print "v0 tp1", distance(v0, tp1)
    print "v1 tp0", distance(v1, tp0)
    print "v1 tp1", distance(v1, tp1)
    if distance(v0, tp0) < dist or distance(v0, tp1) < dist or distance(v1, tp0) < dist or distance(v1, tp1) < dist:
      print "removing line"
      rr, cc, v = line_aa(tmp0[1], tmp0[0], tmp1[1], tmp1[0])
      img[rr,cc] = 0
  return img, max_line_length

    def drawLine(self, centerX, centerY, angle, length, colour):
        """Plot a line from an angle, length and center
        Parameters:
            colour: 0:1 intensity percentage

        Attributes:

        """
        # Calculate the offset of the start and end from the center
        radAngle = math.radians(angle)
        xOffset = length * math.cos(radAngle)
        yOffset = length * math.sin(radAngle)

        # Calculate start and finish of the lines
        startx = math.floor(centerX - xOffset)
        stopx = math.floor(centerX + xOffset)
        starty = math.floor(centerY - yOffset)
        stopy = math.floor(centerY + yOffset)
        # Draw the line onto the array
        rr, cc, val = line_aa(starty, startx, stopy, stopx)
        #rr, cc = line(starty, startx, stopy, stopx)
        self._outputImage[rr, cc] = colour * 255

#.........这里部分代码省略.........

    """
    The third stage aims to produce an encoding that is sensitive to
    local image content while remaining resistant to small changes in
    pose or appearance. The adopted method pools gradient orientation
    information locally in the same way as the SIFT [Lowe 2004]
    feature. The image window is divided into small spatial regions,
    called "cells". For each cell we accumulate a local 1-D histogram
    of gradient or edge orientations over all the pixels in the
    cell. This combined cell-level 1-D histogram forms the basic
    "orientation histogram" representation. Each orientation histogram
    divides the gradient angle range into a fixed number of
    predetermined bins. The gradient magnitudes of the pixels in the
    cell are used to vote into the orientation histogram.
    """

    magnitude = sqrt(gx ** 2 + gy ** 2)
    orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90

    sy, sx = image.shape
    cx, cy = pixels_per_cell
    bx, by = cells_per_block

    n_cellsx = int(np.floor(sx // cx))  # number of cells in x
    n_cellsy = int(np.floor(sy // cy))  # number of cells in y

    # compute orientations integral images
    orientation_histogram = np.zeros((n_cellsy, n_cellsx, orientations))
    for i in range(orientations):
        #create new integral image for this orientation
        # isolate orientations in this range

        temp_ori = np.where(orientation < 180 / orientations * (i + 1),
                            orientation, 0)
        temp_ori = np.where(orientation >= 180 / orientations * i,
                            temp_ori, 0)
        # select magnitudes for those orientations
        cond2 = temp_ori > 0
        temp_mag = np.where(cond2, magnitude, 0)

        orientation_histogram[:,:,i] = uniform_filter(temp_mag, size=(cy, cx))[cy/2::cy, cx/2::cx]


    # now for each cell, compute the histogram
    #orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations))

    radius = min(cx, cy) // 2 - 1
    hog_image = None
    if visualise:
        hog_image = np.zeros((sy, sx), dtype=float)

    if visualise:
        from skimage import draw

        for x in range(n_cellsx):
            for y in range(n_cellsy):
                for o in range(orientations):
                    centre = tuple([y * cy + cy // 2, x * cx + cx // 2])
                    dx = radius * cos(float(o) / orientations * np.pi)
                    dy = radius * sin(float(o) / orientations * np.pi)

                    rr, cc, val = draw.line_aa(centre[0] - int(dx), centre[1] - int(dy), centre[0] + int(dx), centre[1] + int(dy))
                    hog_image[rr, cc] += orientation_histogram[y, x, o]

    """
    The fourth stage computes normalisation, which takes local groups of
    cells and contrast normalises their overall responses before passing
    to next stage. Normalisation introduces better invariance to illumination,
    shadowing, and edge contrast. It is performed by accumulating a measure
    of local histogram "energy" over local groups of cells that we call
    "blocks". The result is used to normalise each cell in the block.
    Typically each individual cell is shared between several blocks, but
    its normalisations are block dependent and thus different. The cell
    thus appears several times in the final output vector with different
    normalisations. This may seem redundant but it improves the performance.
    We refer to the normalised block descriptors as Histogram of Oriented
    Gradient (HOG) descriptors.
    """

    n_blocksx = (n_cellsx - bx) + 1
    n_blocksy = (n_cellsy - by) + 1
    normalised_blocks = np.zeros((n_blocksy, n_blocksx,
                                  by, bx, orientations))

    for x in range(n_blocksx):
        for y in range(n_blocksy):
            block = orientation_histogram[y:y + by, x:x + bx, :]
            eps = 1e-5
            normalised_blocks[y, x, :] = block / sqrt(block.sum() ** 2 + eps)

    """
    The final step collects the HOG descriptors from all blocks of a dense
    overlapping grid of blocks covering the detection window into a combined
    feature vector for use in the window classifier.
    """

    if visualise:
        return normalised_blocks.ravel(), hog_image
    else:
        return normalised_blocks.ravel()

img[rr, cc, :] = (1, 0, 1)
rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=-math.pi / 4.)
img[rr, cc, :] = (0, 0, 1)
rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 2.)
img[rr, cc, :] = (1, 1, 1)

ax1.imshow(img)
ax1.set_title('No anti-aliasing')
ax1.axis('off')


from skimage.draw import line_aa, circle_perimeter_aa


img = np.zeros((100, 100), dtype=np.double)

# anti-aliased line
rr, cc, val = line_aa(12, 12, 20, 50)
img[rr, cc] = val

# anti-aliased circle
rr, cc, val = circle_perimeter_aa(60, 40, 30)
img[rr, cc] = val


ax2.imshow(img, cmap=plt.cm.gray, interpolation='nearest')
ax2.set_title('Anti-aliasing')
ax2.axis('off')

plt.show()

def add_line(arr, p1x, p1y, p2x, p2y):
    from skimage.draw import line_aa
    rr, cc, val = line_aa(p1y, p1x, p2y, p2x)
    arr[rr, cc] = val * 255
    
    return arr

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt

from skimage.transform import hough_line, hough_line_peaks
from skimage.draw import line_aa


if __name__ == '__main__':
    # Construct test image
    image = np.zeros((200, 100))

    # Line 0
    rr, cc, val = line_aa(50, 15, 140, 97)
    image[rr, cc] = val

    # Line 1
    rr, cc, val = line_aa(100, 5, 12, 84)
    image[rr, cc] = np.fmax(image[rr, cc], val)

    # Line 2
    rr, cc, val = line_aa(140, 15, 160, 80)
    image[rr, cc] = np.fmax(image[rr, cc], val)

    # Compute the Hough transform for display
    h, theta, d = hough_line(image)

    fig, ax = plt.subplots(1, 3, figsize=(6, 3.5))

    # Original image
    ax[0].imshow(image, cmap=plt.cm.gray)

def fast_triangle_mesh_drawer(tris, origin, look):
  # shouldn't do anything
  look /= norm(look)

  img = np.zeros((Y, X))
  zimg = np.ones((Y, X))*np.inf

  def project_point(pt):
    # vector from pt to origin
    vv = pt - origin
    v = vv/norm(vv)

    # real projection shit
    vx = npa((v[0], 0, v[2]))
    lx = npa((look[0], 0, look[2]))
    vy = npa((0, v[1], v[2]))
    ly = npa((0, look[1], look[2]))
    def ang(v1, v2):
      v1 /= norm(v1)
      v2 /= norm(v2)
      angl = np.dot(v1, v2)
      crs = np.cross(v1, v2)
      if np.sum(crs) >= 0.0:
        return np.arccos(angl)
      else:
        return -np.arccos(angl)

    x = (ang(vx, lx) / arcrad_per_pixel)
    y = (ang(vy, ly) / arcrad_per_pixel)

    # add z for z-buffering
    # z is the distance of the point from the plane formed by look and origin

    # project v on to look
    z = np.dot(v, look) * norm(vv)

    """
    print " *** "
    print v, K
    print pt, x, y
    """

    return int(round(x + X/2)),int(round(Y/2 - y)), z

  # project the triangles into 2D space
  # does this projection preserve the u and v
  DRAW_WIREFRAME = False
  if DRAW_WIREFRAME:
    lines = []
    for tr in tris:
      p0, p1, p2 = project_point(tr[0]), project_point(tr[1]), project_point(tr[2])
      lines.append((p0[0:2], p1[0:2]))
      lines.append((p1[0:2], p2[0:2]))
      lines.append((p2[0:2], p0[0:2]))

    for pt1, pt2 in lines:
      rr, cc, val = line_aa(pt1[1], pt1[0], pt2[1], pt2[0])

      # filter
      rr[np.logical_or(rr < 0, rr >= Y)] = 0
      cc[np.logical_or(cc < 0, cc >= X)] = 0

      img[rr, cc] = val
  else:
    # z-buffering, keeping the quality low in line with the rest of the program
    polys = []
    for tr in tris:
      xyz = npa(map(project_point, tr))

      # get min and max
      xmin, xmax = np.min(xyz[:, 0]), np.max(xyz[:, 0])
      ymin, ymax = np.min(xyz[:, 1]), np.max(xyz[:, 1])

      # on screen
      xmin = np.clip(xmin, 0, X).astype(np.int)
      xmax = np.clip(xmax, 0, X).astype(np.int)
      ymin = np.clip(ymin, 0, Y).astype(np.int)
      ymax = np.clip(ymax, 0, Y).astype(np.int)

      # triangle in 3 space
      vs1 = xyz[1][0:2] - xyz[0][0:2]
      vs2 = xyz[2][0:2] - xyz[0][0:2]
      vsx = np.cross(vs1, vs2)

      # shade
      shade = random.uniform(0.3, 1.0)

      for x in range(xmin, xmax):
        for y in range(ymin, ymax):
          q = npa([x,y]) - xyz[0][0:2]
          u = np.cross(q, vs2) / vsx
          v = np.cross(vs1, q) / vsx
          if u >= 0 and v >= 0 and u+v <= 1.0:
            pt_xyz = (1-u-v)*xyz[0] + u*xyz[1] + v*xyz[2]
            if pt_xyz[2] < zimg[y,x]:
              zimg[y,x] = pt_xyz[2]
              img[y,x] = shade
      #polys.append((np.mean(xyz[:, 2]), xyz))
     
    """
#.........这里部分代码省略.........