# Create portraits from photos
# GNU GPL 2008


import Image  
import math
import sys
import random
import ImageDraw


def makePortrait(img):
    # This routine implements
    # Sobel edge detection, median smoothing
    # and global thresholding

    # Config parameters       
    maskThresh = 0.85
    maskSize = 3
    # number of smoothing rounds
    numIter = 4
    # magnification for x and y
    xmag = 2
    ymag = 2
    skipSmoothing = False    
    # 
    brush = "brush4"
    harLen = 2
    binaryImage = False
    localThresh = ((maskSize*maskSize) -1)
	
    if img.mode != "RGB":
        img = img.convert("RGB")
    imgdata = img.load()

    imgX = img.size[0] 
    imgY = img.size[1]
    tmpArr1 = map( lambda x: [x]*imgY, [0]*imgX )
    tmpArr2 = map( lambda x: [x]*imgY, [0]*imgX )

    # Sobel kernel
    gx = [[-1, 0, 1], 
          [-2, 0, 2], 
	  [-1, 0, 1]]
    gy = [[ 1, 2, 1], 
	  [ 0, 0, 0], 
	  [-1,-2,-1]]

    count = 0
    mean = 0
    histGlob = [0]*256	    
    hMask = int(maskSize/2)
    for x in xrange(hMask, imgX-hMask):
       for y in xrange(hMask, imgY-hMask):
          pixelGx = pixelGy = 0
	  for i in xrange(-hMask,hMask+1):	    
	     for j in xrange(-hMask,hMask+1):
                pixel = sum(imgdata[x+i,y+j])/3
                pixelGx += gx[i+1][j+1] * pixel
                pixelGy += gy[i+1][j+1] * pixel
          newpixel = math.sqrt(pixelGx * pixelGx + pixelGy * pixelGy)
          newpixel = 255 - int(min(max(newpixel,0),255))
          tmpArr1[x][y] = newpixel
          if skipSmoothing:
             histGlob[pixel] += 1 
             count += 1		  

    inp = tmpArr1
    out = tmpArr2

    outImage1 = Image.new(img.mode, (imgX*xmag,imgY*ymag) , "white")
    outdata1 = outImage1.load()
    draw = ImageDraw.Draw(outImage1)  
    xmax = imgX*xmag
    ymax = imgY*ymag
    xmin = ymin = 0

    if skipSmoothing:
       numIter = 1
       globalThresh = maskThresh * count    	  
       for i in xrange(0, 255):
          mean += histGlob[i]	     
	  #print i,histGlob[i], mean
          if (mean > globalThresh):       
             break	    	    	         
          mean = i
    for g in xrange(0, numIter):
       for x in xrange(maskSize, imgX-maskSize):
          for y in xrange(maskSize, imgY-maskSize):
	     pixel = inp[x][y]	        	    		    
	     if not skipSmoothing:
                hist = [0]*256	    
                mean1 = 0
	        for i in xrange(-maskSize,maskSize+1):	    
	           for j in xrange(-maskSize,maskSize+1):	    
                      if i == 0 and j == 0:
                         continue
                      hist[inp[x+i][y+j]] += 1
                for i in xrange(0, 255):
                   mean1 += hist[i]
                   if (mean1 > localThresh):
                      pixel = i
	    	      break	    	    	    
	     if (x%100 == 0 and y%100 == 0):
	        print g,x,y
	     out[x][y] = pixel   	     
	     if g == numIter-2:             
                histGlob[pixel] += 1 
                count += 1
             if (g == numIter-1) and (pixel < mean):		     	          	       
	          x1 = x*xmag
                  y1 = y*ymag	 
		  if binaryImage:
	             pixel = 0		  
                  outdata1[x1, y1] = (pixel, pixel, pixel)     
		  if brush == "brush1":
                     ang = random.random() * 360
		  if brush == "brush2": 
                     ang = 90
		  if  brush == "brush1" or brush == "brush2":    
		     x2 = min(max(x1 + (harLen * math.cos(ang)),0), xmax-1)
                     y2 = min(max(y1 + (harLen * math.sin(ang)),0), ymax-1)
		     draw.line((x1,y1,x2,y2), fill=(pixel,pixel,pixel))
		  if brush == "brush3":             			 
		     r = harLen
                     draw.ellipse((x1-r,y1-r,x1+r,y1+r), fill=(pixel,pixel,pixel))
		  if brush == "brush4":
                     r = harLen
		     draw.rectangle((x1-r,y1-r,x1+r,y1+r), fill=(pixel,pixel,pixel))
       temp = inp  	    
       inp = out
       out = temp
       if  g == numIter-2:
          globalThresh = maskThresh * count    	  
          for i in xrange(0, 255):
             mean += histGlob[i]	     
	     #print i,histGlob[i], mean
             if (mean > globalThresh):       
                break	    	    	         
          mean = i
          print "global thresh = ", mean, globalThresh, count
    return outImage1



try:
   import psyco
   psyco.full()
except ImportError:
   pass
img = Image.open(sys.argv[1])  
print img.size
outimg = makePortrait(img)  
outimg.save("output.png")