A two-pass paletted pixel-scaling algorithm that uses weighting counting of adjacent colors and a fitness function (for tie-breaking) to create a 2x scale image. This is not an efficient implementation, just a quick-and-dirty proof of concept.

readme_zoom2x.md

The "zoom2x" algorithm

by Andrew G. Crowell

A two-pass paletted pixel-scaling algorithm that uses weighting counting of adjacent colors and a fitness function (for tie-breaking) to create a 2x scale image.

This is not an efficient implementation, just a quick-and-dirty proof of concept. So it is mainly useful for offline rendering right now, but a few optimizations to create less temporary memory and it could be made pretty quick. In particular, the best_sample function will create a dictionary every call, resulting in a lot of garbage. This algorithm could directly work on an indexed image instead and then the weight array be a fixed-length array that is the size of the image color palette (possibly 16 or 256-color or whatever) that's shared between calls and just cleared before use, and then this should result in way fewer allocations. Also somebody could write it in a systems language like C++ or Rust instead of Python -- which would also help a lot, and hopefully wouldn't be too bad to port.

Turns an image like this:

into something like this, 2x as big:

Personally I prefer nearest-neighbour most of the time for scaling game artwork because it gives the cleanest/sharpest look, but I wanted this upscaling algorithm for the purpose of making larger sprites with smaller ones, while maintaining most of the pixel-integrity and original color palette when scaling.

Here's an example of the algorithm in use in pico8 0.16:

These sample images are copyright (c) 2016 Andrew G. Crowell. All rights reserved.

However, I don't mind you using the code for your own projects. I've released the code under the MIT License. I've attached a Python 2.7 version (that uses PIL), and a pico8 0.16 version.

Inspired by the Scale2x algorithm I tried out a new approach. It scales an image to twice the size without interpolation like scale2x, but uses a different way of determining the "best color" for each of the upscaled pixels, which could be thought of a weighted statistical mode of adjacent pixels for each corner.

for a source image of size w x h,
a destination image of size 2w x 2h is created.

for each pixel in the source we sample the pixel p and its neighbours:

[ ul ][ u  ][ ur ]
[ l  ][ p  ][ r  ]
[ dl ][ d  ][ dr ]

A nearest neighbour 2x scaling would output a 2x2 grid of pixels p, like this.
[p][p]
[p][p]

On the other hand, zoom2x will output a 2x2 grid of pixels like this:

[e][f]
[g][h]

e = best_sample(u, l, ul, p)
f = best_sample(u, r, ur, p)
g = best_sample(d, l, dl, p)
h = best_sample(d, r, dr, p)

best_sample(vert, horz, diag, pixel) is a function that is responsible for weighting these pixels.
e, f, g, h essentially are a function of the source pixel, 2 orthogonal adjacent pixels, and a diagonal adjacent pixel.
It uses a combination of a weighted mode, and a fitness function for breaking ties.

The weighting is such that common colors among (horz, vert, pixel) are favored, the source pixel is favored is there is a disagreement, and if a diagonal penalty is provided, diag subtracts from the weight table -- to prevent excess staircasing.

The fitness function for tie-breaking will favor colors around middle brightness, over extremely dark and extremely bright colors because it seemed to give better results (so luminance = (r + g + b) / 3 that is closest to 127).

 I am handwaving like crazy because really I just tried so crap until it worked, but yeah. maybe it's useful to someone. Read the code for more details.
 
 Feel free to make a newer, better, more efficient scaling algorithm if this inspires you and the other ones out there don't fill the void in your life for nicely preserved color/details when upscaling pixels. zzz

zoom2x.lua

function div(n,m)
 return flr(n/m)
end

function unhex(s,n)
 n=n or 2
 local t={}
 for i=1,#s,n do
  add(t,('0x'..sub(s,i,i+n-1))+0)
 end
 return t
end

clrfit=unhex'00010101010102010202010202020201'
samp={}
function bestc(a,b,c,d)
 for i=0,15 do samp[i]=0 end
 samp[a]+=1
 samp[b]+=1
 samp[c]=max(samp[c]-2,0)
 samp[d]+=2
 local best=0
 for i=1,15 do
  if(samp[i]>samp[best] or samp[i]==samp[best] and clrfit[i+1]>clrfit[best+1])best=i
 end
 return best
end

function z2x_(sf,w,h,df,fn)
 function px(x,y)
  if(x<0 or y<0 or x>=w or y>=h)return 0
  return fn(sget(sf%16*8+x,div(sf,16)*8+y))
 end
 for i=0,w-1 do
  for j=0,h-1 do
   local l=px(i-1,j)
   local r=px(i+1,j) 
   local u=px(i,j-1)
   local d=px(i,j+1)
   local p=px(i,j)   
   local dx=df%16*8+i*2
   local dy=div(df,16)*8+j*2
   local e=bestc(u,l,px(i-1,j-1),p)
   local f=bestc(u,r,px(i+1,j-1),p)
   local g=bestc(d,l,px(i-1,j+1),p)
   local h=bestc(d,r,px(i+1,j+1),p)  
   if(e~=0)sset(dx,dy,e)
   if(f~=0)sset(dx+1,dy,f)
   if(g~=0)sset(dx,dy+1,g)
   if(h~=0)sset(dx+1,dy+1,h)
  end
 end
end

-- clear the sprite
function sclr(f,w,h)
 for x=0,w-1 do
  for y=0,h-1 do
   sset(f%16*8+x,div(f,16)*8+y,0)
  end
 end
end

-- scale a sprite in source frame sf of size w,h
-- and store at dest frame df of size w*2,h*2
function z2x(sf,w,h,df)
 sclr(df,w*2,h*2)
 z2x_(sf,w,h,df,function(x)return x==5 and 0 or x end)
 z2x_(sf,w,h,df,function(x)return x~=5 and 0 or x end)
end

for i=0,7 do
 z2x(64+i,8,8,96+i*2)
 z2x(72+i,8,8,128+i*2)
end

zoom2x.py

#!/usr/bin/env python
#
# The "zoom2x algorithm".
#
# by Andrew G. Crowell (@eggboycolor)
#
# Playing around with a method to scale images to 2x size,
# Works slightly better (for me) than scale2x / EPX for my personal use.
#
# This code is released under the MIT License.
#
# Copyright (c) 2016 Andrew G. Crowell
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import PIL.Image

USAGE = '''
source: name of the source file
dest: name of the destination file
sx, sy, sw, sh:
    subregion of the source image (for scaling a single sprite in a spritesheet)
    Defaults to whole image.
diagonal_penalty:
    When non-zero, apply a penalty to adjacent colors that also appear on the
    diagonal adjacent cell. Good for very small sprites (eg. 8x8), bad effects
    on whole-screen scaling.
pico8:
    When non-zero, use a lookup table based on the 16-color pico8 color table
    rather than luminance for masking and fitness. This is more faithful to an
    algorithm I wrote specifically for Pico8, assuming fixed color indices in
    the pico8 palette and using a lookup table rather than calculating luminance.
'''

def best_sample(vert, horz, diag, pixel, diagonal_penalty, fitness):
    weights = {}
    weights[vert] = weights.get(vert, 0) + 1
    weights[horz] = weights.get(horz, 0) + 1
    weights[diag] = max(weights.get(diag, 0) - diagonal_penalty, 0)
    weights[pixel] = weights.get(pixel, 0) + 2
    best = None
    for i in weights.keys():
        if best is None \
        or weights[i] > weights[best] \
        or weights[i] == weights[best] and fitness(i) > fitness(best):
            best = i
    return best if best is not None else (0, 0, 0, 0)

def read_pixel(imagedata, sx, sy, sw, sh, ox, oy, mask):
    if ox < 0 or oy < 0 or ox >= sw or oy >= sh:
        return (0, 0, 0, 0)
    return imagedata[sx + ox, sy + oy] if mask(imagedata[sx + ox, sy + oy]) else (0, 0, 0, 0)

def clear_image(imagedata, x, y, w, h):
    for i in range(w):
        for j in range(h):
            imagedata[x + i, y + j] = (0, 0, 0, 0)

def zoom2x_(sourcedata, sx, sy, sw, sh, destdata, dx, dy, mask, diagonal_penalty, fitness):
    for i in range(sw):
        for j in range(sh):
            l = read_pixel(sourcedata, sx, sy, sw, sh, i - 1, j, mask)
            r = read_pixel(sourcedata, sx, sy, sw, sh, i + 1, j, mask) 
            u = read_pixel(sourcedata, sx, sy, sw, sh, i, j - 1, mask)
            d = read_pixel(sourcedata, sx, sy, sw, sh, i, j + 1, mask)
            p = read_pixel(sourcedata, sx, sy, sw, sh, i, j, mask)
            x = dx + i * 2
            y = dy + j * 2
            e = best_sample(u, l, read_pixel(sourcedata, sx, sy, sw, sh, i - 1, j - 1, mask), p, diagonal_penalty, fitness)
            f = best_sample(u, r, read_pixel(sourcedata, sx, sy, sw, sh, i + 1, j - 1, mask), p, diagonal_penalty, fitness)
            g = best_sample(d, l, read_pixel(sourcedata, sx, sy, sw, sh, i - 1, j + 1, mask), p, diagonal_penalty, fitness)
            h = best_sample(d, r, read_pixel(sourcedata, sx, sy, sw, sh, i + 1, j + 1, mask), p, diagonal_penalty, fitness)  
            if e[3] > 0:
                destdata[x, y] = e
            if f[3] > 0:
                destdata[x + 1, y] = f
            if g[3] > 0:
                destdata[x, y + 1] = g
            if h[3] > 0:
                destdata[x + 1, y + 1] = h

def zoom2x(sourcedata, sx, sy, sw, sh, destdata, dx, dy, mask1, mask2, diagonal_penalty, fitness):
    clear_image(destdata, dx, dy, sw * 2, sh * 2)
    zoom2x_(sourcedata, sx, sy, sw, sh, destdata, dx, dy, mask1, diagonal_penalty, fitness)
    zoom2x_(sourcedata, sx, sy, sw, sh, destdata, dx, dy, mask2, diagonal_penalty, fitness)

def color_luminance(c):
    return (c[0] + c[1] + c[2]) / 3    

def mask_isolate_outline(c):
    return color_luminance(c) < 128

def mask_remove_outline(c):    
    return color_luminance(c) >= 128

def mask_identity(c):    
    return c

def fitness_luminance(c):
    return abs(color_luminance(c) - 128) * (c[3] / 255.0)

def color_distance(c, c2):
    return ((c[0] - c2[0]) * (c[0] - c2[0])
        + (c[1] - c2[1]) * (c[1] - c2[1])
        + (c[2] - c2[2]) * (c[2] - c2[2])) ** 0.5

PICO8_PALETTE = [
    (0x00, 0x00, 0x00),
    (0x1D, 0x2B, 0x53),
    (0x7E, 0x25, 0x53),
    (0x00, 0x87, 0x51),
    (0xAB, 0x52, 0x36),
    (0x5F, 0x57, 0x4F),
    (0xC2, 0xC3, 0xC7),
    (0xFF, 0xF1, 0xE8),
    (0xFF, 0x00, 0x4D),
    (0xFF, 0xA3, 0x00),
    (0xFF, 0xF0, 0x24),
    (0x00, 0xE7, 0x56),
    (0x29, 0xAD, 0xFF),
    (0x83, 0x76, 0x9C),
    (0xFF, 0x77, 0xA8),
    (0xFF, 0xCC, 0xAA),
]
    
PICO8_PALETTE_INDEX = dict((v, i) for i, v in enumerate(PICO8_PALETTE))
PICO8_FITNESS = [0, 1, 1, 1, 1, 1, 2, 1, 2, 2, 1, 2, 2, 2, 2, 1]

def get_pico8_palette_index(c):
    return PICO8_PALETTE_INDEX[min(PICO8_PALETTE, key=lambda i: color_distance(i, c))] if c[3] > 0 else 0

def mask_isolate_outline_pico8(c):
    return get_pico8_palette_index(c) in (1, 2, 5) 

def mask_remove_outline_pico8(c):
    return get_pico8_palette_index(c) not in (1, 2, 5)

def fitness_pico8(c):
    return PICO8_FITNESS[get_pico8_palette_index(c)]

if __name__ == '__main__':
    import sys
    import os.path
    
    if len(sys.argv) < 3:
        exit('Usage: ' + os.path.basename(sys.argv[0]) + ' source dest [sx sy sw sh] [diagonal_penalty] [pico8]\n' + USAGE)
    else:
        source = PIL.Image.open(sys.argv[1])
        source = source.convert('RGBA')
        sourcedata = source.load()
        if len(sys.argv) >= 7:
            x = min(max(0, int(sys.argv[3])), source.size[0] - 1)
            y = min(max(0, int(sys.argv[4])), source.size[1] - 1)
            w = min(max(0, int(sys.argv[5])), source.size[0] - x)
            h = min(max(0, int(sys.argv[6])), source.size[1] - y)
        else:
            x, y, w, h = 0, 0, source.size[0], source.size[1]

        diagonal_penalty = 2 if len(sys.argv) >= 8 and int(sys.argv[7]) != 0 else 0
        pico8 = len(sys.argv) >= 9 and int(sys.argv[8]) != 0

        dest = PIL.Image.new('RGBA', (w * 2, h * 2))
        destdata = dest.load()
        if pico8:
            zoom2x(sourcedata, x, y, w, h, destdata, 0, 0, mask_remove_outline_pico8, mask_isolate_outline_pico8, diagonal_penalty, fitness_pico8)            
        else:
            zoom2x(sourcedata, x, y, w, h, destdata, 0, 0, mask_remove_outline, mask_isolate_outline, diagonal_penalty, fitness_luminance)
        dest.save(open(sys.argv[2], 'wb'))