2

I map spatial data onto a one dimensional interval. First I use a space filling Lebesgue curve (or Z curve) in order to connect my points. This works and I get the following plot if I use gnuplot:

enter image description here

Then I want to transform the Lebesgue curve to a Hilbert curve. But it does not work. This is the output:

enter image description here

So it seems to work at the beginning. But after a while strange things happen and I don't know why.

Here is my code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <limits.h>
#include <stdint.h>

#define DIM 2

// direction of hilbert curve
const unsigned char DirTable[4][4] = 
{{1,2,0,0},{0,1,3,1},{2,0,2,3},{3,3,1,2}};

// Map z-order to hilbert curve
const unsigned char HilbertTable[4][4] =
{{0,3,1,2},{0,1,3,2},{2,3,1,0},{2,1,3,0}};

unsigned int Lebesgue2Hilbert(unsigned int lebesgue){
    unsigned int hilbert = 1;
    int level = 0;
    int dir = 0;
    for(unsigned int tmp=lebesgue; tmp>1; tmp>>=DIM, level++);
    for(; level>0; level--){
        int cell = (lebesgue >> ((level-1)*DIM)) & ((1<<DIM)-1);
        hilbert = (hilbert<<DIM) + HilbertTable[dir][cell];
        dir = DirTable[dir][cell];
    }
    return hilbert;
}

unsigned int Part1By1(unsigned int x)
{
  x &= 0x0000ffff;
  x = (x ^ (x <<  8)) & 0x00ff00ff;
  x = (x ^ (x <<  4)) & 0x0f0f0f0f;
  x = (x ^ (x <<  2)) & 0x33333333;
  x = (x ^ (x <<  1)) & 0x55555555;
  return x;
}

unsigned int EncodeMorton2(unsigned int x, unsigned int y)
{
  return (Part1By1(y) << 1) + Part1By1(x);
}

struct sortKey{
    int idx;
    int key;
};

int compare (const void * a, const void * b)
{
    return ( (*(struct sortKey*)b).key - (*(struct sortKey*)a).key );
}

int uniform_distribution(int rangeLow, int rangeHigh) {
    double myRand = rand()/(1.0 + RAND_MAX); 
    int range = rangeHigh - rangeLow + 1;
    int myRand_scaled = (myRand * range) + rangeLow;
    return myRand_scaled;
}

int main (int argc, char **argv){
    srand(time(NULL));

    int n = 20;
    int N = n*n;

    // Uniform distribution    
    unsigned int pospar[DIM*N];
    int k = 0;
    for(unsigned int i=0; i<n; i++){
        for(unsigned int j=0; j<n; j++){
            pospar[DIM*k] = i;
            pospar[DIM*k+1] = j;
            k++;
        }
    }

    // Lebesgue curve
    unsigned int lkey[N];
    for(int i=0; i<N; i++){
        lkey[i] = EncodeMorton2(pospar[i*DIM+0],pospar[i*DIM+1]);
    //    printf("Lkey: %d\n", lkey[i]);
    }

    // Hilbert curve 
    unsigned int hkey[N];
    for(int i=0; i<N; i++){
        hkey[i] = Lebesgue2Hilbert(lkey[i]);
    }

    struct sortKey sk[N];
    for(int i=0;i<N; i++){
        sk[i].idx = i;
        sk[i].key = hkey[i]; // "lkey[i]" or "hkey[i]"
    }

    qsort (sk, N, sizeof(struct sortKey), compare);

    for (int i=0; i<N; i++){
        printf ("%d %d\n", pospar[DIM*sk[i].idx], pospar[DIM*sk[i].idx+1]);
    }

    return 0;
}

Does anybody see what I have missed?

ad absurdum
  • 19,498
  • 5
  • 37
  • 60
Gilfoyle
  • 3,282
  • 3
  • 47
  • 83
  • Questions seeking debugging help ("why isn't this code working?") must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Questions without a clear problem statement are not useful to other readers. See: How to create a Minimal, Complete, and Verifiable example. – too honest for this site Jan 12 '17 at 20:12
  • 1
    I don't know if it is relevant but in `compare()` supplied to `qsort()` you should not use a one-line subtractive solution. The subtraction can overflow and distort the result. You should make two specific `>` and `<` comparisons with a final `return 0;` – Weather Vane Jan 12 '17 at 20:15
  • @Samuel:Your z order is also wrong. – Micromega Feb 10 '17 at 09:40

0 Answers0