src/trace/filterset.cpp

	filterset.cpp revision a00d232ef8ecb04fb29c1c7ea67f8e14df15b309
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/*
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * Some filters for Potrace in Inkscape
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * Authors:
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *   Bob Jamison <rjamison@titan.com>
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * Copyright (C) 2004 Bob Jamison
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * Released under GNU GPL, read the file 'COPYING' for more information
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#include <cstdio>
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#include <stdlib.h>
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#include "imagemap-gdk.h"
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#include "filterset.h"
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/*#########################################################################
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor### G A U S S I A N  (smoothing)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#########################################################################*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorstatic int gaussMatrix[] =
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     2,  4,  5,  4, 2,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     4,  9, 12,  9, 4,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     5, 12, 15, 12, 5,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     4,  9, 12,  9, 4,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     2,  4,  5,  4, 2
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor};
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorGrayMap *grayMapGaussian(GrayMap *me)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int width  = me->width;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int height = me->height;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstX = 2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastX  = width-3;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstY = 2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastY  = height-3;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *newGm = GrayMapCreate(width, height);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!newGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int y = 0 ; y<height ; y++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        for (int x = 0 ; x<width ; x++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        /* image boundaries */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            if (x<firstX || x>lastX || y<firstY || y>lastY)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                newGm->setPixel(newGm, x, y, me->getPixel(me, x, y));
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                continue;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            /* all other pixels */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            int gaussIndex = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            unsigned long sum = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            for (int i= y-2 ; i<=y+2 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                for (int j= x-2; j<=x+2 ; j++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    int weight = gaussMatrix[gaussIndex++];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sum += me->getPixel(me, j, i) * weight;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            sum /= 159;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        newGm->setPixel(newGm, x, y, sum);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return newGm;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorRgbMap *rgbMapGaussian(RgbMap *me)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int width  = me->width;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int height = me->height;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstX = 2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastX  = width-3;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstY = 2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastY  = height-3;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    RgbMap *newGm = RgbMapCreate(width, height);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!newGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int y = 0 ; y<height ; y++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        for (int x = 0 ; x<width ; x++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        /* image boundaries */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            if (x<firstX || x>lastX || y<firstY || y>lastY)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                newGm->setPixelRGB(newGm, x, y, me->getPixel(me, x, y));
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                continue;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            /* all other pixels */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            int gaussIndex = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            int sumR       = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            int sumG       = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            int sumB       = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            for (int i= y-2 ; i<=y+2 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                for (int j= x-2; j<=x+2 ; j++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    int weight = gaussMatrix[gaussIndex++];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    RGB rgb = me->getPixel(me, j, i);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sumR += weight * (int)rgb.r;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sumG += weight * (int)rgb.g;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sumB += weight * (int)rgb.b;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            RGB rout;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            rout.r = ( sumR / 159 ) & 0xff;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            rout.g = ( sumG / 159 ) & 0xff;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            rout.b = ( sumB / 159 ) & 0xff;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        newGm->setPixelRGB(newGm, x, y, rout);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return newGm;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/*#########################################################################
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor### C A N N Y    E D G E    D E T E C T I O N
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#########################################################################*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorstatic int sobelX[] =
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    -1,  0,  1 ,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    -2,  0,  2 ,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    -1,  0,  1
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor};
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorstatic int sobelY[] =
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     1,  2,  1 ,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor     0,  0,  0 ,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    -1, -2, -1
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor};
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * Perform Sobel convolution on a GrayMap
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorstatic GrayMap *grayMapSobel(GrayMap *gm,
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor               double dLowThreshold, double dHighThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int width  = gm->width;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int height = gm->height;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstX = 1;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastX  = width-2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int firstY = 1;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int lastY  = height-2;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *newGm = GrayMapCreate(width, height);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!newGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int y = 0 ; y<height ; y++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        for (int x = 0 ; x<width ; x++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            unsigned long sum = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        /* image boundaries */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            if (x<firstX || x>lastX || y<firstY || y>lastY)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                sum = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /* ### SOBEL FILTERING #### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                long sumX = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                long sumY = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                int sobelIndex = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                for (int i= y-1 ; i<=y+1 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    for (int j= x-1; j<=x+1 ; j++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        sumX += gm->getPixel(gm, j, i) *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             sobelX[sobelIndex++];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                sobelIndex = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                for (int i= y-1 ; i<=y+1 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    for (int j= x-1; j<=x+1 ; j++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        sumY += gm->getPixel(gm, j, i) *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             sobelY[sobelIndex++];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*###  GET VALUE ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                sum = abs(sumX) + abs(sumY);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (sum > 765)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sum = 765;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#if 0
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*###  GET ORIENTATION (slow, pedantic way) ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                double orient = 0.0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (sumX==0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    if (sumY==0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        orient = 0.0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    else if (sumY<0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        sumY = -sumY;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        orient = 90.0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        orient = 90.0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    orient = 57.295779515 * atan2( ((double)sumY),((double)sumX) );
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    if (orient < 0.0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        orient += 180.0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*###  GET EDGE DIRECTION ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                int edgeDirection = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (orient < 22.5)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    edgeDirection = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            else if (orient < 67.5)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    edgeDirection = 45;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            else if (orient < 112.5)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    edgeDirection = 90;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            else if (orient < 157.5)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    edgeDirection = 135;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*###  GET EDGE DIRECTION (fast way) ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                int edgeDirection = 0; /*x,y=0*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (sumX==0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    if (sumY!=0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        edgeDirection = 90;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   /*long slope = sumY*1024/sumX;*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   long slope = (sumY << 10)/sumX;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   if (slope > 2472 || slope< -2472)  /*tan(67.5)*1024*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                       edgeDirection = 90;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   else if (slope > 414) /*tan(22.5)*1024*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                       edgeDirection = 45;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   else if (slope < -414) /*-tan(22.5)*1024*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                       edgeDirection = 135;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                   }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#endif
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /* printf("%ld %ld %f %d\n", sumX, sumY, orient, edgeDirection); */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*### Get two adjacent pixels in edge direction ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                unsigned long leftPixel;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                unsigned long rightPixel;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (edgeDirection == 0)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    leftPixel  = gm->getPixel(gm, x-1, y);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    rightPixel = gm->getPixel(gm, x+1, y);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else if (edgeDirection == 45)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    leftPixel  = gm->getPixel(gm, x-1, y+1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    rightPixel = gm->getPixel(gm, x+1, y-1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else if (edgeDirection == 90)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    leftPixel  = gm->getPixel(gm, x, y-1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    rightPixel = gm->getPixel(gm, x, y+1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else /*135 */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    leftPixel  = gm->getPixel(gm, x-1, y-1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    rightPixel = gm->getPixel(gm, x+1, y+1);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*### Compare current value to adjacent pixels ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                /*### if less that either, suppress it ### */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                if (sum < leftPixel || sum < rightPixel)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    sum = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    unsigned long highThreshold =
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                          (unsigned long)(dHighThreshold * 765.0);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    unsigned long lowThreshold =
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                          (unsigned long)(dLowThreshold * 765.0);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    if (sum >= highThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        sum = 765; /* EDGE.  3*255 this needs to be settable */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    else if (sum < lowThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        sum = 0; /* NONEDGE */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        if ( gm->getPixel(gm, x-1, y-1)> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x  , y-1)> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x+1, y-1)> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x-1, y  )> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x+1, y  )> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x-1, y+1)> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x  , y+1)> highThreshold ||
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                             gm->getPixel(gm, x+1, y+1)> highThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                            sum = 765; /* EDGE fix me too */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                            sum = 0; /* NONEDGE */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                }/* else */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            if (sum==0) /* invert light & dark */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                sum = 765;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            else
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                sum = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            newGm->setPixel(newGm, x, y, sum);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        }/* for (x) */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    }/* for (y) */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return newGm;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorGrayMap *
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorgrayMapCanny(GrayMap *gm, double lowThreshold, double highThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!gm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *gaussGm = grayMapGaussian(gm);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!gaussGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    /*gaussGm->writePPM(gaussGm, "gauss.ppm");*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *cannyGm = grayMapSobel(gaussGm, lowThreshold, highThreshold);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!cannyGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    /*cannyGm->writePPM(cannyGm, "canny.ppm");*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    gaussGm->destroy(gaussGm);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return cannyGm;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorGdkPixbuf *
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorgdkCanny(GdkPixbuf *img, double lowThreshold, double highThreshold)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!img)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *grayMap = gdkPixbufToGrayMap(img);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!grayMap)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    /*grayMap->writePPM(grayMap, "gbefore.ppm");*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GrayMap *cannyGm = grayMapCanny(grayMap,lowThreshold, highThreshold);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    grayMap->destroy(grayMap);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!cannyGm)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    /*grayMap->writePPM(grayMap, "gafter.ppm");*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    GdkPixbuf *newImg = grayMapToGdkPixbuf(cannyGm);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return newImg;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/*#########################################################################
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor### Q U A N T I Z A T I O N
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor#########################################################################*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylortypedef struct OctreeNode_def OctreeNode;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * The basic octree node
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorstruct OctreeNode_def
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned long r;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned long g;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned long b;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned int  index;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned long nrPixels;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    unsigned int  nrChildren;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    OctreeNode *parent;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    OctreeNode *children[8];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor};
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor * create an octree node, and initialize it
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill TaylorOctreeNode *octreeNodeCreate()
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    OctreeNode *node = (OctreeNode *)malloc(sizeof(OctreeNode));
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!node)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->r             = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->g             = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->b             = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->index         = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->nrPixels      = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->nrChildren    = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->parent        = NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int i=0 ; i<8 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->children[i] = NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return node;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *  delete an octree node and its children
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorvoid octreeNodeDelete(OctreeNode *node)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!node)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int i=0 ; i<8 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        octreeNodeDelete(node->children[i]);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    free(node);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *  delete the children of an octree node
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorvoid octreeNodeDeleteChildren(OctreeNode *node)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    if (!node)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        return;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    node->nrChildren = 0;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int i=0 ; i<8 ; i++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        octreeNodeDelete(node->children[i]);
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->children[i]=NULL;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor/**
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *  insert an RGB value into an octree node according to its
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor *  high-order rgb vector bits
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylorint octreeNodeInsert(OctreeNode *root, RGB rgb, int bitsPerSample)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor{
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    OctreeNode *node = root;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int newColor     = FALSE;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int r            = rgb.r;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int g            = rgb.g;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int b            = rgb.b;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    int shift = 7;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    for (int bit=0 ; bit<bitsPerSample ; bit++)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        /* update values of all nodes from the root to the leaf */
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->r += r;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->g += g;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->b += b;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node->nrPixels++;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        int index = (((r >> shift) & 1) << 2) |
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                (((g >> shift) & 1) << 1) |
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor                    (((b >> shift) & 1)     ) ;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        OctreeNode *child = node->children[index];
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        if (!child)
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            {
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            child                 = octreeNodeCreate();
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            node->children[index] = child;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            child->parent         = node;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            node->nrChildren++;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        newColor              = TRUE;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor            }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        node = child; /*next level*/
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        shift--;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor        }
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor    return newColor;
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor}
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor
9e39c5ba00a55fa05777cc94b148296af305e135Bill Taylor

/**
 *  find an exact match for an RGB value, at the given bits
 *  per sample.  if not found, return -1
 */
int octreeNodeFind(OctreeNode *root, RGB rgb, int bitsPerSample)
{
    OctreeNode *node = root;
    int r            = rgb.r;
    int g            = rgb.g;
    int b            = rgb.b;

    int shift = 7;
    for (int bit=0 ; bit<bitsPerSample ; bit++)
        {
        int index = (((r >> shift) & 1) << 2) |
                (((g >> shift) & 1) << 1) |
                    (((b >> shift) & 1)     ) ;

        OctreeNode *child = node->children[index];
        if (!child)
            return -1;
        node = child; /*next level*/
        shift--;
        }
    printf("error.  this should not happen\n");
    return -1;
}


static void spaces(int nr)
{
    for (int i=0; i<nr ; i++)
        printf(" ");
}

/**
 *  pretty-print an octree node and its children
 */
void octreeNodePrint(OctreeNode *node, int indent)
{
    spaces(indent); printf("####Node###\n");
    spaces(indent); printf("r :%lu\n", node->r);
    spaces(indent); printf("g :%lu\n", node->g);
    spaces(indent); printf("b :%lu\n", node->b);
    spaces(indent); printf("i :%d\n", node->index);
    for (unsigned int i=0; i<8; i++)
        {
        OctreeNode *child = node->children[i];
        if (!child)
            continue;
        spaces(indent); printf("   child %d :", i);
        octreeNodePrint(child, indent+4);
        }
}

/**
 *  Count how many nodes in this (sub)tree
 */
static int octreeNodeCount(OctreeNode *node)
{
    int count = 1;
    for (unsigned int i=0; i<8; i++)
        {
        OctreeNode *child = node->children[i];
        if (!child)
            continue;
        count += octreeNodeCount(child);
        }
    return count;
}


/**
 * Count all of the leaf nodes in the octree
 */
static void octreeLeafArray(OctreeNode *node, OctreeNode **array, int arraySize, int *len)
{
    if (!node)
        return;
    if (node->nrChildren == 0 && *len < arraySize)
        {
        array[*len] = node;
        *len = *len + 1;
        }
    for (int i=0 ; i<8 ; i++)
        octreeLeafArray(node->children[i], array, arraySize, len);
}


/**
 *  Count all of the leaf nodes in the octree
 */
static int octreeLeafCount(OctreeNode *node)
{
    if (!node)
        return 0;
    if (node->nrChildren == 0)
        return 1;
    int leaves = 0;
    for (int i=0 ; i<8 ; i++)
        leaves += octreeLeafCount(node->children[i]);
    return leaves;
}

/**
 * Mark all of the leaf nodes in the octree with an index nr
 */
static void octreeLeafIndex(OctreeNode *node, int *index)
{
    if (!node)
        return;
    if (node->nrChildren == 0)
        {
        node->index = *index;
        *index = *index + 1;
        return;
        }
    for (int i=0 ; i<8 ; i++)
        octreeLeafIndex(node->children[i], index);
}


/**
 * Find a node that has children, and that also
 * has the lowest pixel count
 */
static void octreefindLowestLeaf(OctreeNode *node, OctreeNode **lowestLeaf)
{
    if (!node)
        return;
    if (node->nrChildren == 0)
        {
        if (node->nrPixels < (*lowestLeaf)->nrPixels)
            *lowestLeaf = node;
        return;
        }

    for (int i=0 ; i<8 ; i++)
        octreefindLowestLeaf(node->children[i], lowestLeaf);
}


/**
 * reduce the leaves on an octree to a given number
 */
int octreePrune(OctreeNode *root, int nrColors)
{
    int leafCount = octreeLeafCount(root);

    while (leafCount > nrColors)
        {
        OctreeNode *lowestLeaf = root;
        octreefindLowestLeaf(root, &lowestLeaf);

        if (!lowestLeaf)
            break; //should never happen

       if (lowestLeaf==root)
            {
            printf("found no leaves\n");
            continue;
            }

        OctreeNode *parent = lowestLeaf->parent;
        if (!parent)
            continue;

        for (int i=0 ; i<8 ; i++)
            {
            OctreeNode *child = parent->children[i];
            if (child == lowestLeaf)
                {
                parent->nrChildren--;
                octreeNodeDelete(child);
                parent->children[i] = NULL;
                break;
                }
            }
        /*printf("leafCount:%d lowPixels:%lu\n",
               leafCount, lowestLeaf->nrPixels);*/
        leafCount = octreeLeafCount(root);
        }
    int index = 0;
    octreeLeafIndex(root, &index);

    //printf("leafCount:%d\n", leafCount);
    //octreeNodePrint(root, 0);

    return leafCount;
}


/**
 *
 */
OctreeNode *octreeBuild(RgbMap *rgbMap, int bitsPerSample, int nrColors)
{
    OctreeNode *root = octreeNodeCreate();
    if (!root)
        return NULL;
    for (int y=0 ; y<rgbMap->height ; y++)
        {
        for (int x=0 ; x<rgbMap->width ; x++)
            {
            RGB rgb = rgbMap->getPixel(rgbMap, x, y);
            octreeNodeInsert(root, rgb, bitsPerSample);
            }
        }

    if (octreePrune(root, nrColors)<0)
        {
        octreeNodeDelete(root);
        return NULL;
        }

    /* octreeNodePrint(root, 0); */

    return root;
}


/* Compare two rgb's for brightness, for the qsort() below */
static int compRGB(const void *a, const void *b)
{
    RGB *ra = (RGB *)a;
    RGB *rb = (RGB *)b;
    int ba = ra->r + ra->g + ra->b;
    int bb = rb->r + rb->g + rb->b;
    int comp = ba - bb;
    return comp;
}

/**
 *
 */
RGB *makeRGBPalette(OctreeNode *root, int nrColors)
{

    OctreeNode **palette = (OctreeNode **)malloc(nrColors * sizeof(OctreeNode *));
    if (!palette)
        {
        return NULL;
        }
    int len = 0;
    octreeLeafArray(root, palette, nrColors, &len);

    RGB *rgbpal = (RGB *)malloc(len * sizeof(RGB));
    if (!rgbpal)
        {
        free(palette);
        return NULL;
        }

    for (int i=0; i<len ; i++)
        {
        OctreeNode *node = palette[i];
        RGB rgb;
        if (node->nrPixels == 0)
            {
            continue;
            }
        rgb.r = (unsigned char)( (node->r / node->nrPixels) & 0xff);
        rgb.g = (unsigned char)( (node->g / node->nrPixels) & 0xff);
        rgb.b = (unsigned char)( (node->b / node->nrPixels) & 0xff);
        int index = node->index;
        //printf("Pal: %d %d %d %d\n", rgb.r, rgb.g, rgb.b, index);
        rgbpal[index]=rgb;
        }

    free(palette);

    /* sort by brightness, to avoid high contrasts */
    qsort((void *)rgbpal, len, sizeof(RGB), compRGB);

    return rgbpal;
}


/**
 *  Return the closest color in the palette to the request
 */
RGB lookupQuantizedRGB(RGB *rgbpalette, int paletteSize, RGB candidate, int *closestIndex)
{
    /* slow method */
    unsigned long closestMatch = 10000000;
    RGB closestRGB = { 0 , 0, 0 };
    *closestIndex = 0;
    for (int i=0 ; i<paletteSize ; i++)
        {
        RGB entry = rgbpalette[i];
        unsigned long dr    = candidate.r - entry.r;
        unsigned long dg    = candidate.g - entry.g;
        unsigned long db    = candidate.b - entry.b;
        unsigned long match = dr * dr + dg * dg + db * db;
        if (match < closestMatch)
            {
            *closestIndex = i;
            closestMatch  = match;
            closestRGB    = entry;
            }
        }

    return closestRGB;
}


/**
 * Quantize an RGB image to a reduced number of colors.  bitsPerSample
 * is usually 3 - 5 out of 8 to conserve cpu and memory
 */
IndexedMap *rgbMapQuantize(RgbMap *rgbMap, int bitsPerSample, int nrColors)
{
    if (!rgbMap)
        return NULL;

    OctreeNode *otree = octreeBuild(rgbMap, bitsPerSample, nrColors);
    if (!otree)
        {
        return NULL;
        }

    /*Make sure we don't request more colors than actually exist*/
    int nodeCount = octreeLeafCount(otree);
    if (nodeCount < nrColors)
        nrColors = nodeCount;

    RGB *rgbpal = makeRGBPalette(otree, nrColors);
    if (!rgbpal)
        {
        octreeNodeDelete(otree);
        return NULL;
        }

    /*We have our original and palette. Make the new one*/
    IndexedMap *newMap = IndexedMapCreate(rgbMap->width, rgbMap->height);
    if (!newMap)
        {
        free(rgbpal);
        octreeNodeDelete(otree);
        return NULL;
        }

    /* fill in the color lookup table */
    for (int i=0 ; i< nrColors ; i++)
        {
        newMap->clut[i] = rgbpal[i];
        }
    newMap->nrColors = nrColors;

    for (int y=0 ; y<rgbMap->height ; y++)
        {
        for (int x=0 ; x<rgbMap->width ; x++)
            {
            RGB rgb = rgbMap->getPixel(rgbMap, x, y);
            //int indexNr = octreeNodeFind(otree, rgb, bitsPerSample);
            //printf("i:%d\n", indexNr);
            //RGB quantRgb = rgbpal[indexNr];
            int closestIndex;
            RGB quantRgb = lookupQuantizedRGB(rgbpal, nrColors, rgb, &closestIndex);
            newMap->setPixel(newMap, x, y, (unsigned int)closestIndex);
            }
        }

    free(rgbpal);
    octreeNodeDelete(otree);

    return newMap;
}


/**
 *  Experimental.  Work on this later
 */
GrayMap *quantizeBand(RgbMap *rgbMap, int nrColors)
{

    int bitsPerSample = 4;

    RgbMap *gaussMap = rgbMapGaussian(rgbMap);
    //gaussMap->writePPM(gaussMap, "rgbgauss.ppm");

    IndexedMap *qMap = rgbMapQuantize(gaussMap, bitsPerSample, nrColors);
    //qMap->writePPM(qMap, "rgbquant.ppm");
    gaussMap->destroy(gaussMap);

    GrayMap *gm = GrayMapCreate(rgbMap->width, rgbMap->height);

    /* RGB is quantized.  There should now be a small set of (R+G+B) */
    for (int y=0 ; y<qMap->height ; y++)
        {
        for (int x=0 ; x<qMap->width ; x++)
            {
            RGB rgb = qMap->getPixelValue(qMap, x, y);
            int sum = rgb.r + rgb.g + rgb.b;
            if (sum & 1)
                sum = 765;
            else
                sum = 0;
            /*printf("%d %d %d : %d\n", rgb.r, rgb.g, rgb.b, index);*/
            gm->setPixel(gm, x, y, sum);
            }
        }

    qMap->destroy(qMap);

    return gm;
}


/*#########################################################################
### E N D    O F    F I L E
#########################################################################*/