图像渲染

#include "stdafx.h"
#include <SDL.h>
#include <SDL_image.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <float.h>
#include <limits.h>


#ifdef WIN32
#pragma comment(lib, "SDL.lib")
#pragma comment(lib, "SDLmain.lib")
#pragma comment(lib, "SDL_image.lib")
#endif


#define BRUTE       0x0
#define GAUSSIAN    0x1
#define AVG         0x2
#define ANN         0x4


#define YIQ         0x0
#define RGB         0x1


#define DIFF        0x0
#define PSRC        0x1


/* the images we load (or save, in case of g_B) */
static SDL_Surface * g_a = NULL;
static SDL_Surface * g_A = NULL;
static SDL_Surface * g_b = NULL;
static SDL_Surface * g_B = NULL;


/* our own pixel abstraction */
typedef struct 
{
    int x;
    int y;
} pixel;


/* our own point abstraction (this is really a feature vector and
 * probably should be an array, but since I know I'll be only dealing
 * with RGB (or YIQ) I only need to worry about 3 values) 
 */
typedef struct
{
    int x; /* the r or y */
    int y; /* the g or i */
    int z; /* the b or q */
} point;


/* our pyramids */
typedef struct 
{
    point *** points; /* triple pointer, of: level, row, column */
    int * h; /* array of height values per level */
    int * w; /* array of width values per level */
    pixel ** s;
} pyramid;


static pyramid g_pa;
static pyramid g_pb;
static pyramid g_pA;
static pyramid g_pB;


/* some globals */
static int g_seed = 42; /* seed for randomizing g_B's lowest level if wanted */
static int g_neigh_width = 5; /* width of the neighborhood */
static int g_pyramid_levels = 4; /* the number of pyramid levels, 0 indexed */
static int g_pyramid_type = BRUTE; /* the method used to generate higher pyramid 
                               levels */
static int g_rand = 0; /* whether we want to randomize the lowest level of g_B */
static int g_is_graphical = 1; /* whether the program runs in graphical mode */
static int g_color_mode = RGB; /* the color mode we want to run in (RGB requires 
                           3 times as many computations) */
static int g_passes = 1; /* number of passes to do at each pyramid level */
static int g_method_type = BRUTE; /* method to compute neighborhoods */
static int g_source = DIFF; /* method to calculate new pixel */
static double g_coherence = 0.0; /* for coherence checking */


void usage(char * execname);
void handle_opts(int argc, char ** argv);
void draw_images(SDL_Surface * screen);
SDL_Color get_pixel(SDL_Surface * surf, int x, int y);
void set_pixel(SDL_Surface * surf, int x, int y, SDL_Color color);
void draw_progress(SDL_Surface * screen, int pos);
void draw_bar(SDL_Surface * screen, int width);


void do_brute_analogy(SDL_Surface * screen);
pixel best_brute_match(pixel q, int l, int pass);
int do_dist(point a, point b, int weight);
int do_brute_dist(pixel p, pixel q, int l, int pass);
void do_diff(int l, pixel p, pixel q);


/* since the purpose of YIQ is speedup, we write a whole set of
 * essentially duplicate functions so we don't have to hit a boatload
 * of if statements everywhere
 */
void do_brute_analogy_YIQ(SDL_Surface * screen);
pixel best_brute_match_YIQ(pixel q, int l, int pass);
int do_dist_YIQ(point a, point b, int weight);
int do_brute_dist_YIQ(pixel p, pixel q, int l, int pass);
void do_diff_YIQ(int l, pixel p, pixel q);


void create_pyramids(int levels, int type, SDL_Surface * screen);
point get_point(pyramid p, int level, int x, int y);
void pyramid_init_brute(void);
void pyramid_init_avg(void);
void pyramid_init_gaussian(void);
point do_avg(pyramid p, int l, int x, int y);
point do_gauss(pyramid p, int l, int x, int y);


SDL_Color p_get_color(pyramid p, int level, int x, int y);
void draw_pyramid_level(SDL_Surface * screen, pyramid p, int level, int xoff, int yoff);
point to_YIQ(point color);
point to_RGB(point color);
void save_image(void);


int main(int argc, char ** argv)
{
    SDL_Surface * screen;


    fprintf(stdout, "Handling options & initializing...\n");
    handle_opts(argc, argv);
    srand((unsigned int)g_seed);
    fprintf(stdout, "Initializing SDL...\n");
        
    if (SDL_Init(SDL_INIT_VIDEO) != 0)
        perror("SDL Init failed");
    (void)atexit(SDL_Quit);


    if (g_is_graphical==1)
    {
        screen = SDL_SetVideoMode(800, 640, 16, SDL_DOUBLEBUF);


        if(screen == NULL)
            perror("SDL Set Video Mode failed"), exit(EXIT_FAILURE);
    }
    
    /* draw our images */
    fprintf(stdout, "Initializing pyramids\n");
    create_pyramids(g_pyramid_levels, g_pyramid_type, screen);
    fprintf(stdout, "Finished making pyramids\n");


    if (g_is_graphical==1)
    {
        fprintf(stdout, "Drawing initial pyramids\n");
        draw_images(screen);
        SDL_Flip(screen);
        fprintf(stdout, "Done drawing initial pyramids\n");
    }
    
    /* TODO: ANN? */
    
    switch (g_method_type)
    {
    case BRUTE:
        if (g_color_mode == RGB)
            do_brute_analogy(screen);
        else if (g_color_mode == YIQ)
            do_brute_analogy_YIQ(screen);
        break;
    case ANN:
        /* TODO: implement ANN */
        fprintf(stdout, "Sorry ANN currently not supported\n");
        break;
    }


    /* convert lowest level of g_pB to our image and save it */
    save_image();
    
    /* TODO: massive memory cleanup, for now let OS handle it */
    SDL_FreeSurface(g_a);
    SDL_FreeSurface(g_A);
    SDL_FreeSurface(g_b);
    SDL_FreeSurface(g_B);
    SDL_FreeSurface(screen);


    return EXIT_SUCCESS;
}


void
save_image()
{
    int x, y;
    SDL_Color col;


    fprintf(stdout, "Saving image as output.bmp\n");
    for (y=0; y<g_B->h; y++)
    {
        for (x=0; x<g_B->w; x++)
        {
            col = p_get_color(g_pB, 0, x, y);
            set_pixel(g_B, x, y, col);
        }
    }
    fprintf(stdout, "Done...\n");
    SDL_SaveBMP(g_B, "output.bmp");
}
/** 
 * Draws the initial images 
 * @param   screen  The screen we're drawing to.
 */
void
draw_images(SDL_Surface * screen)
{ /* {{{ */
    int l;
    int xoff, yoff;


    /* this won't handle the new YIQ colors :( 
    SDL_Rect src, dst;
    
    src.x = 0;
    src.y = 0;
    src.w = g_a->w;
    src.h = g_a->h;
    dst=src;
    SDL_BlitSurface(g_a, &src, screen, &dst);
   
    src.x = 0;
    src.y = 0;
    src.w = g_A->w;
    src.h = g_A->h;
    dst=src;
    dst.y = g_a->h + 10;
    SDL_BlitSurface(g_A, &src, screen, &dst);
    
    src.x = 0;
    src.y = 0;
    src.w = g_b->w;
    src.h = g_b->h;
    dst=src;
    dst.y = g_a->h+10+g_A->h+10;
    SDL_BlitSurface(g_b, &src, screen, &dst);
    
    src.x = 0;
    src.y = 0;
    src.w = g_B->w;
    src.h = g_B->h;
    dst=src;
    dst.y = g_a->h+10+g_A->h+10+g_b->h+10;


    SDL_BlitSurface(g_B, &src, screen, &dst);
    */
    xoff = 0;
    for (l = 0, yoff=0; l < g_pyramid_levels; l++)
    {
        draw_pyramid_level(screen, g_pa, l, xoff, yoff);
        yoff+=g_pa.h[l]+10;
    }
    xoff += g_pa.w[0] + 10;
    for (l = 0, yoff=0; l < g_pyramid_levels; l++)
    {
        draw_pyramid_level(screen, g_pA, l, xoff, yoff);
        yoff+=g_pA.h[l]+10;
    }
    xoff += g_pA.w[0] + 10;
    for (l = 0, yoff=0; l < g_pyramid_levels; l++)
    {
        draw_pyramid_level(screen, g_pb, l, xoff, yoff);
        yoff+=g_pb.h[l]+10;
    }
    xoff += g_pb.w[0] + 10;
    for (l = 0, yoff=0; l < g_pyramid_levels; l++)
    {
        draw_pyramid_level(screen, g_pB, l, xoff, yoff);
        yoff+=g_pB.h[l]+10;
    }
     
    SDL_Flip(screen);
    return;
} /* }}} */


/**
 * Handles our input options.
 *
 * @param   argc    the number of arguments
 * @param   argv    the actual arguments
 */
void 
handle_opts(int argc, char ** argv)
{ /* {{{ */
    int a;
    int x, y;
    SDL_Color col;
    
    while( (a=getopt(argc, argv, "Ghra:A:b:B:c:k:l:m:n:p:s:S:t:")) != EOF)
    {
        switch(a)
        {
        case 'G':
            g_is_graphical = 0;
            break;
        case 'h':
            usage(argv[0]);
            exit(EXIT_SUCCESS);
            /* NOT REACHED */
            break;
        case 'r':
            g_rand = 1;
            break;
        case 'a':
            fprintf(stdout, "Loading `%s' as image `a'\n", optarg);
            g_a = IMG_Load(optarg);
            break;
        case 'A':
            fprintf(stdout, "Loading `%s' as image `A'\n", optarg);
            g_A = IMG_Load(optarg);
            break;
        case 'b':
            fprintf(stdout, "Loading `%s' as image `b'\n", optarg);
            g_b = IMG_Load(optarg);
            g_B = IMG_Load(optarg);
            col.r = col.b = col.g = 0;
            for (y=0; y < g_B->h; y++)
            {
                for (x=0; x < g_B->w; x++)
                {
                    set_pixel( g_B, x, y, col);
                }
            }
            break;
        case 'B':
            g_B = IMG_Load(optarg);
            fprintf(stdout, "Loading `%s' as image `B'\n", optarg);
            g_rand = 0;
            break;
        case 'c':
            if (strncmp(optarg, "RGB", 3) == 0) g_color_mode = RGB;
            if (strncmp(optarg, "YIQ", 3) == 0) g_color_mode = YIQ;
            break;
        case 'k':
            g_coherence = atof(optarg);
            break;
        case 'l':
            g_pyramid_levels = atoi(optarg);
            break;
        case 'm':
            if (strncmp(optarg, "BRUTE", 5) == 0) g_method_type = BRUTE;
            if (strncmp(optarg, "ANN", 2) == 0) g_method_type = ANN;
            break;
        case 'n':
            g_neigh_width = atoi(optarg);
            break;
        case 'p':
            g_passes = atoi(optarg);
            break;
        case 's':
            g_seed = atoi(optarg);
            break;
        case 'S':
            if (strncmp(optarg, "PSRC", 4) == 0) g_source = PSRC;
            if (strncmp(optarg, "DIFF", 4) == 0) g_source = DIFF;
            break;
        case 't':
            if (strncmp(optarg, "BRUTE", 5) == 0) g_pyramid_type = BRUTE;
            if (strncmp(optarg, "AVG", 3) == 0) g_pyramid_type = AVG;
            if (strncmp(optarg, "GAUSSIAN", 8) == 0) g_pyramid_type = GAUSSIAN;
            break;
        case '?':
            /* 
             * getopt handles error reporting by default 
             */
            exit(EXIT_FAILURE);
            /* NOT REACHED */
        }
    }
    argc-=optind;
    argv+=optind;


    /* check for extra parameters */
    if (argc > 0)
    {
        fprintf(stderr, "ERROR: too many arguments\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_a == NULL)
    {
        fprintf(stderr, "ERROR: image `a' not specified\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_A == NULL)
    {
        fprintf(stderr, "ERROR: image `A' not specified\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_b == NULL)
    {
        fprintf(stderr, "ERROR: image `b' not specified\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_B == NULL)
    {
        fprintf(stderr, "ERROR: image `B' (noise) not specified\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_neigh_width <= 0)
    {
        fprintf(stderr, "ERROR: neighborhood width must be positive!\n");
        exit(EXIT_FAILURE);
        /* NOT REACHED */
    }
    if (g_passes <= 0)
    {
        fprintf(stderr, "ERROR: passes must be positive\n");
        exit(EXIT_FAILURE);
    }
    /* do randomizing if desired */
    if (g_rand != 0)
    {
        for (y=0; y < g_B->h; y++)
        {
            for (x=0; x < g_B->w; x++)
            {
                col.r = (Uint8) (rand() % 256);
                col.b = (Uint8) (rand() % 256);
                col.g = (Uint8) (rand() % 256);
                set_pixel( g_B, x, y, col);
            }
        }
    }


    /* print what we're running with */
    if (g_is_graphical == 0)
        fprintf(stdout, "Running in non-gui mode\n");
    else
        fprintf(stdout, "Running in gui-mode\n");
    
    fprintf(stdout, "Using neighborhood width of %d\n", g_neigh_width);
    fprintf(stdout, "Pyramid Levels are %d\n", g_pyramid_levels);
    switch(g_color_mode)
    {
        case RGB:
            fprintf(stdout, "Running in RGB\n");
            break;
        case YIQ: 
            fprintf(stdout, "Running in YIQ\n");
            break;
    }
    fprintf(stdout, "Seed is %d\n", g_seed);
    fprintf(stdout, "Number of passes per pyramid level %d\n", g_passes);
    switch(g_pyramid_type)
    {
        case BRUTE:
            fprintf(stdout, "Pyramids constructed using BRUTE method\n");
            break;
        case AVG: 
            fprintf(stdout, "Pyramids constructed using AVG method\n");
            break;
        case GAUSSIAN: 
            fprintf(stdout, "Pyramids constructed using GAUSSIAN method\n");
            break;
    }
    if (g_rand==1) fprintf(stdout, "Pyramid for B randomly initializes\n");
    if (g_coherence>0.0) fprintf(stdout, "Doing coherence with k=%f\n",
            g_coherence);
    switch(g_source)
    {
        case DIFF:
            fprintf(stdout, "Matches used to compute difference\n");
            break;
        case PSRC: 
            fprintf(stdout, "Matches directly applied from A (a2) image\n");
            break;
    }


    return;
} /* }}} */


/**
 * Prints out a usage message
 * @param   execname    the name of the executable
 */
void
usage(char * execname)
{ /* {{{ */
    /* To be ISO compliant strings can't be longer than 509 chars */
    fprintf(stdout, 
            "USAGE: %s [options]\n",
            execname);
    fprintf(stdout,
            "[-G] run in non-gui mode\n"
            "[-h] displays this help menu\n"
            "[-P] if using DIFF, later passes than 1 drop multiplier\n"
            "[-a filename] specifies the a image (required)\n"
            "[-A filename] specifies the A image, i.e., changed a (required)\n"
            "[-b filename] specifies the b image (required)\n" 
            "[-B filename] specifies the noise to init the B image as (optional)\n"
            "[-c color] the color mode to run in (RGB, YIQ)\n");
    fprintf(stdout,
            "[-m method] method to search with (BRUTE, ...)\n"
            "[-n neighborhood width] allows to set the neigborhood width (must be odd)\n"
            "[-p passes] the number of passes to perform at each pyramid level\n"
            "[-r] init B with random noise from rand()\n"
            "[-s seed] allows to specify random seed used for the -r option\n\n"
            "[-S source] allows to specify how to treat best matches (DIFF, PSRC)\n\n"
            "[-t type] the way pyramids are constructed (BRUTE, AVG, GAUSSIAN)\n"
           );
} /* }}} */


/**
 * Sets a pixel in the surface (converts from YIQ if needed)
 * @param   surf    the surface we're writing to
 * @param   x       the x coordinate of the pixel to change
 * @param   y       the y coordinate of the pixel to change
 * @param   color   the color to change it to
 */
void 
set_pixel(SDL_Surface * surf, int x, int y, SDL_Color color)
{ /* {{{ */
    Uint32 col;
    char * pos;
    
    col = SDL_MapRGB(surf->format, color.r, color.g, color.b);
    pos = (char *) surf->pixels;
    pos += (surf->pitch * y);
    pos += (surf->format->BytesPerPixel * x);
    memcpy(pos, &col, surf->format->BytesPerPixel);
} /* }}} */


/**
 * Gets a pixel at a particular location (converts to YIQ if needed)
 * @param   surf    the surface we're getting the pixel from
 * @param   x       the x coordinate of the pixel to get
 * @param   y       the y coordinate of the pixel to get
 */
SDL_Color
get_pixel(SDL_Surface * surf, int x, int y)
{ /* {{{ */
    SDL_Color color;
    Uint32 col;
    char * pos;
    
    pos = (char *)surf->pixels;
    pos += (surf->pitch * y);
    pos += (surf->format->BytesPerPixel * x);
    
    memcpy(&col, pos, surf->format->BytesPerPixel);
    SDL_GetRGB(col, surf->format, &color.r, &color.g, &color.b);
    return color;
} /* }}} */


void
do_brute_analogy(SDL_Surface * screen)
{
    int l, x, y;
    int pass;
    pixel p, q;
    SDL_Event event;
    SDL_keysym keysym;
    int xoff, yoff;
    SDL_Color col;


    xoff = g_a->w + 10 + g_A->w + 10 + g_b->w + 10;


    for (l=g_pyramid_levels-1; l>=0; l--)
    {
        yoff = 0;
        for (x=0; x<l; x++)
            yoff+=g_pB.h[x]+10;
        pass = 0;
        for (y=0; y < g_pB.h[l]; y++)
        {
            for (x=0; x < g_pB.w[l]; x++)
            {
                q.x=x;  q.y=y;
                p = best_brute_match(q, l, pass);
                g_pB.points[l][x][y] = g_pA.points[l][p.x][p.y]; 
                if (g_source == DIFF)
                {
                    do_diff(l, p, q);
                }
                col = p_get_color(g_pB, l, x, y);
                if(g_is_graphical==1) set_pixel(screen, xoff + x, yoff+y, col);
                while (SDL_PollEvent(&event) != 0)
                {
                    switch (event.type)
                    {
                        case SDL_KEYDOWN:
                            keysym = event.key.keysym;


                            if (keysym.sym == SDLK_q)
                            {
                                fprintf(stdout, "Exiting\n");
                                save_image();
                                exit(EXIT_SUCCESS);
                            }
                            break;
                    }
                }
            }
            if(g_is_graphical==1)SDL_Flip(screen);
            else fprintf(stdout, "level: %d\tpass: %d\trow: %d of %d\n",
                    l, pass, y, g_pB.h[l]);
        }
        pass++;
        while (pass < g_passes)
        {
            for (y=0; y < g_pB.h[l]; y++)
            {
                for (x=0; x < g_pB.w[l]; x++)
                {
                    q.x=x;  q.y=y;
                    p = best_brute_match(q, l, pass);
                    g_pB.points[l][x][y] = g_pA.points[l][p.x][p.y]; 
                    if (g_source == DIFF)
                    {
                        do_diff(l, p, q);
                    }
                    col = p_get_color(g_pB, l, x, y);
                    if(g_is_graphical==1)set_pixel(screen, xoff + x, yoff+y, col);


                    while (SDL_PollEvent(&event) != 0)
                    {
                        switch (event.type)
                        {
                            case SDL_KEYDOWN:
                                keysym = event.key.keysym;


                                if (keysym.sym == SDLK_q)
                                {
                                    fprintf(stdout, "Exiting\n");
                                    save_image();
                                    exit(EXIT_SUCCESS);
                                }
                                break;
                        }
                    }
                }
                if(g_is_graphical==1)SDL_Flip(screen);
                else fprintf(stdout, "level: %d\tpass: %d\trow: %d of %d\n",
                        l, pass, y, g_pB.h[l]);
            }
            pass++;
        }
    }
    return;
}


void
do_diff(int l, pixel p, pixel q)
{
    g_pB.points[l][q.x][q.y].x = g_pb.points[l][q.x][q.y].x + 
        g_pA.points[l][p.x][p.y].x - g_pa.points[l][p.x][p.y].x;
    g_pB.points[l][q.x][q.y].y = g_pb.points[l][q.x][q.y].y + 
        g_pA.points[l][p.x][p.y].y - g_pa.points[l][p.x][p.y].y;
    g_pB.points[l][q.x][q.y].z = g_pb.points[l][q.x][q.y].z + 
        g_pA.points[l][p.x][p.y].z - g_pa.points[l][p.x][p.y].z;


    /* normalize */
    if (g_pB.points[l][q.x][q.y].x > 255) g_pB.points[l][q.x][q.y].x=255;
    else if(g_pB.points[l][q.x][q.y].x < 0) g_pB.points[l][q.x][q.y].x=0;
    if (g_pB.points[l][q.x][q.y].y > 255) g_pB.points[l][q.x][q.y].y=255;
    else if(g_pB.points[l][q.x][q.y].y < 0) g_pB.points[l][q.x][q.y].y=0;
    if (g_pB.points[l][q.x][q.y].z > 255) g_pB.points[l][q.x][q.y].z=255;
    else if(g_pB.points[l][q.x][q.y].z < 0) g_pB.points[l][q.x][q.y].z=0;
}


pixel
best_brute_match(pixel q, int l, int pass)
{
    pixel p; 
    pixel match; 
    int x, y;
    int best_dist = INT_MAX;
    int dist;
    /*
    int w = g_neigh_width / 2;
    pixel coh_match;
    int coh_dist;
    int best_coh_dist = INT_MAX;
    */
    /* go through all the pixels in example images (a, A) and find the
     * one that minimizes distance
     */
    match.x = 0;
    match.y = 0;
    for (x=0; x < g_pA.w[l]; x++)
    {
        for (y=0; y < g_pA.h[l]; y++)
        {
            p.x=x; p.y=y;
            dist = do_brute_dist(p, q, l, pass);
            if (dist < best_dist)
            {
                match = p;
                best_dist = dist;
            }
        }
    }
    /* TODO: do coherence */
    /*for (y = -w; y < w; y++)
    {
        for (x = -w; x < w; x++)
        {
            p.x=q.x + w; p.y=q.y + y;
            
            coh_dist = do_brute_dist(p, q, l, pass);
            if (coh_dist < best_coh_dist)
            {
                coh_match = p;
                best_coh_dist = dist;
            }
        }
    }*/
    return match;
}


int
do_brute_dist(pixel p, pixel q, int l, int pass)
{
    point a, A;
    point b, B;


    int dist = 0;
    int w = g_neigh_width / 2;
    int x, y;
    
    for (y = -w; y < w; y++)
    {
        for (x = -w; x < w; x++)
        {
            /* dist += get_ab_dist(p.x+x, p.y+y); */
            a = get_point(g_pa, l, p.x + x, p.y + y);
            b = get_point(g_pb, l, q.x + x, q.y + y);


            dist += do_dist(a, b, w);
            if (y > 0 && pass==0) continue;
            if (y == 0 && x >= 0 && pass==0) continue;
            
            A = get_point(g_pA, l, p.x + x, p.y + y);
            B = get_point(g_pB, l, q.x + x, q.y + y);
            
            dist += do_dist(A, B, w);
        }
    }
    /* handle one level up (if it exists) */
    if (l + 1 < g_pyramid_levels)
    {
        w/=2;
        
        for (y = -w; y < w; y++)
        {
            for (x = -w; x < w; x++)
            {
                /* dist += get_ab_dist(p.x+x, p.y+y); */
                a = get_point(g_pa, l+1, p.x/2 + x, p.y/2 + y);
                b = get_point(g_pb, l+1, q.x/2 + x, q.y/2 + y);


                dist += do_dist(a, b, w);


                A = get_point(g_pA, l+1, p.x/2 + x, p.y/2 + y);
                B = get_point(g_pB, l+1, q.x/2 + x, q.y/2 + y);


                dist += do_dist(A, B, w);
            }
        }
    }
    return dist;
}


int 
do_dist(point a, point b, int weight)
{
    return 
            ((a.x - b.x) * (a.x - b.x) +
            (a.y - b.y) * (a.y - b.y) +
            (a.z - b.z) * (a.z - b.z)) >> abs(weight);
}


void
do_brute_analogy_YIQ(SDL_Surface * screen)
{
    int l, x, y;
    int pass;
    pixel p, q;
    SDL_Event event;
    SDL_keysym keysym;
    int xoff, yoff;
    SDL_Color col;


    xoff = g_a->w + 10 + g_A->w + 10 + g_b->w + 10;


    for (l=g_pyramid_levels-1; l>=0; l--)
    {
        yoff = 0;
        for (x=0; x<l; x++)
            yoff+=g_pB.h[x]+10;
        pass = 0;
        for (y=0; y < g_pB.h[l]; y++)
        {
            for (x=0; x < g_pB.w[l]; x++)
            {
                q.x=x;  q.y=y;
                p = best_brute_match_YIQ(q, l, pass);
                g_pB.points[l][x][y] = g_pA.points[l][p.x][p.y]; 
                g_pB.points[l][x][y].y = g_pb.points[l][x][y].y;
                g_pB.points[l][x][y].z = g_pb.points[l][x][y].z;
                if (g_source == DIFF)
                {
                    do_diff_YIQ(l, p, q);
                }
                col = p_get_color(g_pB, l, x, y);
                set_pixel(screen, xoff + x, yoff+y, col);
                while (SDL_PollEvent(&event) != 0)
                {
                    switch (event.type)
                    {
                        case SDL_KEYDOWN:
                            keysym = event.key.keysym;


                            if (keysym.sym == SDLK_q)
                            {
                                fprintf(stdout, "Exiting\n");
                                save_image();
                                exit(EXIT_SUCCESS);
                            }
                            break;
                    }
                }
            }
            SDL_Flip(screen);
        }
        pass++;
        while (pass < g_passes)
        {
            for (y=0; y < g_pB.h[l]; y++)
            {
                for (x=0; x < g_pB.w[l]; x++)
                {
                    q.x=x;  q.y=y;
                    p = best_brute_match_YIQ(q, l, pass);
                    g_pB.points[l][x][y] = g_pA.points[l][p.x][p.y]; 
                    g_pB.points[l][x][y].y = g_pb.points[l][x][y].y;
                    g_pB.points[l][x][y].z = g_pb.points[l][x][y].z;
                    if (g_source == DIFF)
                    {
                        do_diff_YIQ(l, p, q);
                    }
                    col = p_get_color(g_pB, l, x, y);
                    set_pixel(screen, xoff + x, yoff+y, col);


                    while (SDL_PollEvent(&event) != 0)
                    {
                        switch (event.type)
                        {
                            case SDL_KEYDOWN:
                                keysym = event.key.keysym;


                                if (keysym.sym == SDLK_q)
                                {
                                    fprintf(stdout, "Exiting\n");
                                    save_image();
                                    exit(EXIT_SUCCESS);
                                }
                                break;
                        }
                    }
                }
                SDL_Flip(screen);
            }
            pass++;
        }
    }
    return;
}


void
do_diff_YIQ(int l, pixel p, pixel q)
{
    g_pB.points[l][q.x][q.y].x = g_pb.points[l][q.x][q.y].x + 
        g_pA.points[l][p.x][p.y].x - g_pa.points[l][p.x][p.y].x;
    g_pB.points[l][q.x][q.y].y = g_pb.points[l][q.x][q.y].y;
    g_pB.points[l][q.x][q.y].z = g_pb.points[l][q.x][q.y].z; 
    /*
    g_pB.points[l][q.x][q.y].y = g_pb.points[l][q.x][q.y].y + 
        g_pA.points[l][p.x][p.y].y - g_pa.points[l][p.x][p.y].y;
    g_pB.points[l][q.x][q.y].z = g_pb.points[l][q.x][q.y].z + 
        g_pA.points[l][p.x][p.y].z - g_pa.points[l][p.x][p.y].z;
    */


    /* normalize */
    /*
    if (g_pB.points[l][q.x][q.y].x > 255) g_pB.points[l][q.x][q.y].x=255;
    else if(g_pB.points[l][q.x][q.y].x < 0) g_pB.points[l][q.x][q.y].x=0;
    if (g_pB.points[l][q.x][q.y].y > 255) g_pB.points[l][q.x][q.y].y=255;
    else if(g_pB.points[l][q.x][q.y].y < 0) g_pB.points[l][q.x][q.y].y=0;
    if (g_pB.points[l][q.x][q.y].z > 255) g_pB.points[l][q.x][q.y].z=255;
    else if(g_pB.points[l][q.x][q.y].z < 0) g_pB.points[l][q.x][q.y].z=0;
    */
}


pixel
best_brute_match_YIQ(pixel q, int l, int pass)
{
    pixel p; 
    pixel match; 
    int x, y;
    int best_dist = INT_MAX;
    int dist;
    /*
    int w = g_neigh_width / 2;
    pixel coh_match;
    int coh_dist;
    int best_coh_dist = INT_MAX;
    */
    
    /* go through all the pixels in example images (a, A) and find the
     * one that minimizes distance
     */
    match.x = 0;
    match.y = 0;
    for (x=0; x < g_pA.w[l]; x++)
    {
        for (y=0; y < g_pA.h[l]; y++)
        {
            p.x=x; p.y=y;
            dist = do_brute_dist_YIQ(p, q, l, pass);
            if (dist < best_dist)
            {
                match = p;
                best_dist = dist;
            }
        }
    }
    /* do coherence */
    /*
    for (y = -w; y < w; y++)
    {
        for (x = -w; x < w; x++)
        {
            p.x=q.x + w; p.y=q.y + y;
            
            coh_dist = do_brute_dist_YIQ(p, q, l, pass);
            if (coh_dist < best_coh_dist)
            {
                coh_match = p;
                best_coh_dist = dist;
            }
        }
    }*/
    return match;
}


int
do_brute_dist_YIQ(pixel p, pixel q, int l, int pass)
{
    point a, A;
    point b, B;


    int dist = 0;
    int w = g_neigh_width / 2;
    int x, y;
    
    for (y = -w; y < w; y++)
    {
        for (x = -w; x < w; x++)
        {
            /* dist += get_ab_dist(p.x+x, p.y+y); */
            a = get_point(g_pa, l, p.x + x, p.y + y);
            b = get_point(g_pb, l, q.x + x, q.y + y);


            dist += do_dist_YIQ(a, b, w);
            if (y > 0 && pass==0) continue;
            if (y == 0 && x >= 0 && pass==0) continue;
            
            A = get_point(g_pA, l, p.x + x, p.y + y);
            B = get_point(g_pB, l, q.x + x, q.y + y);
            
            dist += do_dist_YIQ(A, B, w);
        }
    }
    /* handle one level up (if it exists) */
    if (l + 1 < g_pyramid_levels)
    {
        w/=2;
        
        for (y = -w; y < w; y++)
        {
            for (x = -w; x < w; x++)
            {
                /* dist += get_ab_dist(p.x+x, p.y+y); */
                a = get_point(g_pa, l+1, p.x/2 + x, p.y/2 + y);
                b = get_point(g_pb, l+1, q.x/2 + x, q.y/2 + y);


                dist += do_dist_YIQ(a, b, w);


                A = get_point(g_pA, l+1, p.x/2 + x, p.y/2 + y);
                B = get_point(g_pB, l+1, q.x/2 + x, q.y/2 + y);


                dist += do_dist_YIQ(A, B, w);
            }
        }
    }
    return dist;
}


int 
do_dist_YIQ(point a, point b, int weight)
{
    return ((a.x - b.x) * (a.x - b.x))>> abs(weight);
}


/** Create pyramids according to specified method */
void
create_pyramids(int levels, int type, SDL_Surface * screen)
{ /* {{{ */
    int x, y;
    SDL_Color col;
    
    /* allocate the heights and widths of our levels */
    if ((g_pa.h = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pb.h = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pA.h = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pB.h = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pa.w = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pb.w = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pA.w = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    if ((g_pB.w = (int *)malloc(sizeof(int) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    
    /* init height and width info for all levels */
    for (x=0; x<levels; x++)
    {
        g_pa.h[x] = g_a->h >> x;
        g_pa.w[x] = g_a->w >> x;
        g_pb.h[x] = g_b->h >> x;
        g_pb.w[x] = g_b->w >> x;
        g_pA.h[x] = g_A->h >> x;
        g_pA.w[x] = g_A->w >> x;
        g_pB.h[x] = g_B->h >> x;
        g_pB.w[x] = g_B->w >> x;
    }
    
    /* allocate pointers to the appropriate levels */
    if ((g_pa.points = (point ***)malloc(sizeof(point **) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    for (x=0; x < levels; x++)
    {
        /* allocate the appropriate number of rows at each level */
        if((g_pa.points[x] = (point **)
                    malloc(sizeof(point *) * g_pa.w[0])) == NULL)
            perror("malloc"), exit(EXIT_FAILURE);
        /* allocate the appropriate number of columns at each level */
        for (y=0; y < g_pa.w[x]; y++)
            if((g_pa.points[x][y] = (point *)
                        malloc(sizeof(point) * g_pa.h[0])) == NULL)
                perror("malloc"), exit(EXIT_FAILURE);
    }
    if ((g_pb.points = (point ***)malloc(sizeof(point **) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    for (x=0; x < levels; x++)
    {
        if((g_pb.points[x] = (point **)
                    malloc(sizeof(point *) * g_pb.w[x])) == NULL)
            perror("malloc"), exit(EXIT_FAILURE);
        for (y=0; y < g_pb.w[x]; y++)
            if((g_pb.points[x][y] = (point *)
                        malloc(sizeof(point) * g_pb.h[x])) == NULL)
                perror("malloc"), exit(EXIT_FAILURE);
    }
    if ((g_pA.points = (point ***)malloc(sizeof(point **) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    for (x=0; x < levels; x++)
    {
        if((g_pA.points[x] = (point **)
                    malloc(sizeof(point *) * g_pA.w[x])) == NULL)
            perror("malloc"), exit(EXIT_FAILURE);
        for (y=0; y < g_pA.w[x]; y++)
            if((g_pA.points[x][y] = (point *)
                        malloc(sizeof(point) * g_pA.h[x])) == NULL)
                perror("malloc"), exit(EXIT_FAILURE);
    }
    if ((g_pB.points = (point ***)malloc(sizeof(point **) * levels)) == NULL)
        perror("malloc"), exit(EXIT_FAILURE);
    for (x=0; x < levels; x++)
    {
        if((g_pB.points[x] = (point **)
                    malloc(sizeof(point *) * g_pB.w[x])) == NULL)
            perror("malloc"), exit(EXIT_FAILURE);
        for (y=0; y < g_pB.w[x]; y++)
            if((g_pB.points[x][y] = (point *)
                        malloc(sizeof(point) * g_pB.h[x])) == NULL)
                perror("malloc"), exit(EXIT_FAILURE);
    }
        
    /* populate highest resolution (lowest level) of pyramids */
    for (y=0; y < g_pa.h[0]; y++)
    {
        for (x=0; x < g_pa.w[0]; x++)
        {
            col = get_pixel(g_a, x, y);
            g_pa.points[0][x][y].x = col.r;
            g_pa.points[0][x][y].y = col.g;
            g_pa.points[0][x][y].z = col.b;
            if (g_color_mode == YIQ)
                g_pa.points[0][x][y] = to_YIQ(g_pa.points[0][x][y]);
            /*
            set_pixel(screen, x, y, col);
            col = p_get_color(g_pa, 0, x, y);
            col.r = (Uint8) g_pa.points[0][x][y].x;
            col.g = (Uint8) g_pa.points[0][x][y].y;
            col.b = (Uint8) g_pa.points[0][x][y].z;
            set_pixel(screen, x+200, y, col);
        }
            SDL_Flip(screen);
            SDL_Delay(1);
            for (b=0; b < g_pa.h[0]; b++)
            {
                for (a=0; a <g_pa.w[0]; a++)
                {
                    col.r = (Uint8) g_pa.points[0][a][b].x;
                    col.g = (Uint8) g_pa.points[0][a][b].y;
                    col.b = (Uint8) g_pa.points[0][a][b].z;
                    set_pixel(screen, a+400, b, col);


                }
            }
            SDL_Flip(screen);*/
        }
    }
    for (y=0; y < g_b->h; y++)
    {
        for (x=0; x < g_b->w; x++)
        {
            col = get_pixel(g_b, x, y);
            g_pb.points[0][x][y].x = col.r;
            g_pb.points[0][x][y].y = col.g;
            g_pb.points[0][x][y].z = col.b;
            if (g_color_mode == YIQ)
                g_pb.points[0][x][y] = to_YIQ(g_pb.points[0][x][y]);
        }
    }
    for (y=0; y < g_A->h; y++)
    {
        for (x=0; x < g_A->w; x++)
        {
            col = get_pixel(g_A, x, y);
            g_pA.points[0][x][y].x = col.r;
            g_pA.points[0][x][y].y = col.g;
            g_pA.points[0][x][y].z = col.b;
            if (g_color_mode == YIQ)
                g_pA.points[0][x][y] = to_YIQ(g_pA.points[0][x][y]);
        }
    }
    for (y=0; y < g_B->h; y++)
    {
        for (x=0; x < g_B->w; x++)
        {
            col = get_pixel(g_B, x, y);
            g_pB.points[0][x][y].x = col.r;
            g_pB.points[0][x][y].y = col.g;
            g_pB.points[0][x][y].z = col.b;
            if (g_color_mode == YIQ)
                g_pB.points[0][x][y] = to_YIQ(g_pB.points[0][x][y]);
        }
    }


    /* create higher level pyramids */
    switch(type)
    {
    case BRUTE:
        pyramid_init_brute();
        break;
    case GAUSSIAN:
        pyramid_init_gaussian();
        break;
    case AVG:
        pyramid_init_avg();
        break;
    }
} /* }}} */


point
to_YIQ(point color)
{
    point tmp;


    tmp.x = (int)((.299 * color.x) + (.587 * color.y) + (.114 * color.z));
    tmp.y = (int)((.596 * color.x) - (.274 * color.y) - (.322 * color.z));
    tmp.z = (int)((.212 * color.x) - (.523 * color.y) + (.311 * color.z));
    if (tmp.x < 0) tmp.x = 0;
    else if(tmp.x > 255) tmp.x = 255;
    if (tmp.y < 0) tmp.y = 0;
    else if(tmp.y > 255) tmp.y = 255;
    if (tmp.z < 0) tmp.z = 0;
    else if(tmp.z > 255) tmp.z = 255;


    return tmp;
}


point
to_RGB(point color)
{
    point tmp;


    tmp.x = (int) (color.x + ( .956 * color.y) + ( .621 * color.z));
    tmp.y = (int) (color.x - ( .272 * color.y) - ( .647 * color.z));
    tmp.z = (int) (color.x - (1.105 * color.y) + (1.702 * color.z));


    return tmp;
}


/** initializes the upper pyramids by simply taking the corresponding
 * pixel at the upper level.  Dumb, and loses info.
 */
void
pyramid_init_brute(void)
{ /* {{{ */
    int l, x, y;
    int levels = g_pyramid_levels;


    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pa.h[l]; y++)
        {
            for (x=0; x < g_pa.w[l]; x++)
            {
                g_pa.points[l][x][y]=g_pa.points[l-1][x<<1][y<<1];
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pb.h[l]; y++)
        {
            for (x=0; x < g_pb.w[l]; x++)
            {
                g_pb.points[l][x][y]=g_pb.points[l-1][x<<1][y<<1];
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pA.h[l]; y++)
        {
            for (x=0; x < g_pA.w[l]; x++)
            {
                g_pA.points[l][x][y]=g_pA.points[l-1][x<<1][y<<1];
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pB.h[l]; y++)
        {
            for (x=0; x < g_pB.w[l]; x++)
            {
                g_pB.points[l][x][y]=g_pB.points[l-1][x<<1][y<<1];
            }
        }
    }
} /* }}} */


/** initializes the upper pyramids by simply taking the corresponding
 * average of the "cross" of pixels at the upper level.  
 */
void
pyramid_init_avg(void)
{ /* {{{ */
    int l, x, y;
    int levels = g_pyramid_levels;
    
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pa.h[l]; y++)
        {
            for (x=0; x < g_pa.w[l]; x++)
            {
                g_pa.points[l][x][y] = do_avg(g_pa, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pb.h[l]; y++)
        {
            for (x=0; x < g_pb.w[l]; x++)
            {
                g_pb.points[l][x][y] = do_avg(g_pb, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pA.h[l]; y++)
        {
            for (x=0; x < g_pA.w[l]; x++)
            {
                g_pA.points[l][x][y] = do_avg(g_pA, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pB.h[l]; y++)
        {
            for (x=0; x < g_pB.w[l]; x++)
            {
                g_pB.points[l][x][y] = do_avg(g_pB, l, x, y);
            }
        }
    }
} /* }}} */


point
do_avg(pyramid p, int l, int x, int y)
{ /* {{{ */
    point pts[9];
    point pt = {0, 0, 0};
    int lcv;
    
    pts[0] = get_point(p, l-1, x<<1     , y<<1);
    pts[1] = get_point(p, l-1, (x<<1)+1 , y<<1);
    pts[2] = get_point(p, l-1, x<<1     , (y<<1)+1);
    pts[3] = get_point(p, l-1, (x<<1)-1 , y<<1);
    pts[4] = get_point(p, l-1, x<<1     , (y<<1)-1);
    pts[5] = get_point(p, l-1, (x<<1)-1 , (y<<1)-1);
    pts[6] = get_point(p, l-1, (x<<1)+1 , (y<<1)-1);
    pts[7] = get_point(p, l-1, (x<<1)+1 , (y<<1)+1);
    pts[8] = get_point(p, l-1, (x<<1)-1 , (y<<1)+1);


    for(lcv = 0; lcv < 9; lcv++)
    {
        pt.x+=(int)(pts[lcv].x * .1/.9);
        pt.y+=(int)(pts[lcv].y * .1/.9);
        pt.z+=(int)(pts[lcv].z * .1/.9);
    }
    return pt;
} /* }}} */


/** initializes the upper pyramids by simply taking the corresponding
 * gaussian distribution of the "cross" of pixels at the upper level.  
 */
void
pyramid_init_gaussian(void)
{ /* {{{ */
    int l, x, y;
    int levels = g_pyramid_levels;
    
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pa.h[l]; y++)
        {
            for (x=0; x < g_pa.w[l]; x++)
            {
                g_pa.points[l][x][y] = do_gauss(g_pa, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pb.h[l]; y++)
        {
            for (x=0; x < g_pb.w[l]; x++)
            {
                g_pb.points[l][x][y] = do_gauss(g_pb, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pA.h[l]; y++)
        {
            for (x=0; x < g_pA.w[l]; x++)
            {
                g_pA.points[l][x][y] = do_gauss(g_pA, l, x, y);
            }
        }
    }
    for (l=1; l<levels; l++)
    {
        for (y=0; y < g_pB.h[l]; y++)
        {
            for (x=0; x < g_pB.w[l]; x++)
            {
                g_pB.points[l][x][y] = do_gauss(g_pB, l, x, y);
            }
        }
    }
} /* }}} */


point
do_gauss(pyramid p, int l, int x, int y)
{ /* {{{ */
    point pts[9];
    point pt = {0, 0, 0};
    int lcv;
    
    pts[0] = get_point(p, l-1, x<<1     , y<<1);
    pts[1] = get_point(p, l-1, (x<<1)+1 , y<<1);
    pts[2] = get_point(p, l-1, x<<1     , (y<<1)+1);
    pts[3] = get_point(p, l-1, (x<<1)-1 , y<<1);
    pts[4] = get_point(p, l-1, x<<1     , (y<<1)-1);
    pts[5] = get_point(p, l-1, (x<<1)-1 , (y<<1)-1);
    pts[6] = get_point(p, l-1, (x<<1)+1 , (y<<1)-1);
    pts[7] = get_point(p, l-1, (x<<1)+1 , (y<<1)+1);
    pts[8] = get_point(p, l-1, (x<<1)-1 , (y<<1)+1);


    for(lcv = 1; lcv < 9; lcv++)
    {
        pt.x+=(int)(pts[lcv].x *.08);
        pt.y+=(int)(pts[lcv].y *.08);
        pt.z+=(int)(pts[lcv].z *.08);
    }
    pt.x+=(int)(pts[0].x *.36);
    pt.y+=(int)(pts[0].y *.36);
    pt.z+=(int)(pts[0].z *.36);
    return pt;
} /* }}} */


/** gets a point at a particular level in our pyramid, correctly
 * wrapping around for values that are too big or too small.
 */
point
get_point(pyramid p, int level, int x, int y)
{
    if (x < 0) 
    {
        for(; x < 0; x+=p.w[level]);
    }
    else if (x >= p.w[level]) x%= p.w[level];


    if (y < 0) 
    {
        for(; y < 0; y+=p.h[level]);
    }
    else if (y >= p.h[level]) y%= p.h[level];


    return p.points[level][x][y];
}


SDL_Color
p_get_color(pyramid p, int level, int x, int y)
{
    point pt;
    SDL_Color col;
    
    pt = get_point(p, level, x, y);
    if (g_color_mode == YIQ)
        pt = to_RGB(pt);
    col.r = (Uint8) (pt.x & 0x00FF);
    col.g = (Uint8) (pt.y & 0x00FF);
    col.b = (Uint8) (pt.z & 0x00FF);
    return col;
}


void
draw_pyramid_level(SDL_Surface * screen, pyramid p, int level, int xoff, int yoff)
{
    int x, y; 
    SDL_Color col;
    for (y=0; y<p.h[level]; y++)
    {
        for (x=0; x<p.w[level]; x++)
        {
            col = p_get_color(p, level, x, y);
            set_pixel(screen, x+xoff, y+yoff, col);
        }
    }
}