0

我试图运行用 C 编写的希尔伯特曲线代码，我在这里找到了 http://www.tddft.org/svn/octopus/trunk/src/grid/hilbert.c

代码运行，但我从输出中得到的结果不正确。我制作了一个简单的驱动程序例程，它从命令行获取 3 个值作为参数，并将它们传递给希尔伯特曲线编码、解码例程。更准确地说，我无法解码回原始坐标 (x,y,z)。

我的问题之一是了解变量nbits在做什么。我假设它是编码希尔伯特值的大小。为了检查这一点，我尝试修改其中一个函数的原始定义

void InttoTranspose(const int dim, const long long int h, int * x)

至

void InttoTranspose(const int dim, const long long int h, int * x, int* size)

我分配给*size位计数变量 ifbit 的位置。无论如何，这一切都没有奏效。因此，我想寻求您的帮助。

修改后的代码在这里：

#include<stdio.h>

//#include <config.h>

#include <assert.h>

/* This code is based on the implementation of the Hilbert curves
   presented in:

   J. Skilling, Programming the Hilbert curve, AIP Conf. Proc. 707, 381 (2004); http://dx.doi.org/10.1063/1.1751381

*/

*/ The int* size was included later in an attempt to find the proper size  /*
/* void InttoTranspose(const int dim, const long long int h, int * x)*/

void InttoTranspose(const int dim, const long long int h, int * x, int* size){
  /* the code uses some funny way of storing the bits */

  int idir, ibit, ifbit;

  for(idir = 0; idir < dim; idir++) x[idir] = 0;

  ifbit = 0;
  for(ibit = 0; ibit < 21; ibit++){
    for(idir = dim - 1; idir >= 0; idir--){
      x[idir] += (((h>>ifbit)&1)<<ibit);
      ifbit++;
    }
  }
 *size=ifbit; // I think that this should be nbits
}


void TransposetoAxes(int* X, int b, int n ){ /* position, #bits, dimension */

  int N = 2 << (b-1), P, Q, t;
  int i;

  /* Gray decode by H ^ (H/2) */
  t = X[n-1] >> 1;
  for(i = n - 1; i > 0; i--) X[i] ^= X[i-1];
  X[0] ^= t;

  /* Undo excess work */
  for( Q = 2; Q != N; Q <<= 1 ) {
    P = Q - 1;
    for( i = n-1; i >= 0 ; i-- ){
      if( X[i] & Q ) {
    X[0] ^= P; /* invert */
      } else{ 
    t = (X[0]^X[i]) & P; X[0] ^= t; X[i] ^= t; /* exchange */
      }
    }
  }

}

//void FC_FUNC_(hilbert_index_to_point, HILBERT_INDEX_TO_POINT)(const int * dim, const int * nbits, const long long int * index, int * point){

 // InttoTranspose(*dim, *index, point);
 // TransposetoAxes(point, *nbits, *dim);

//}


int main(int argc,char* argv[]){

  long long int hilbert;
  int i=0,x[2],m;

  while(argc--){
    if(m=atoi(*argv++)) x[i++]=m, printf("--> %5d %5d\n",m,argc);
  }


  printf("x= %5d y= %5d  z= %5d \n",x[0],x[1],x[2]);

  InttoTranspose(3, hilbert, x,&m);
  printf("hilbert encoded --> %llu size -->%d \n",hilbert,m);
  TransposetoAxes(x,m, 3 );
  printf("x= %5d y= %5d  z= %5d \n",x[0],x[1],x[2]);

  return 0;
}
4

0 回答 0