File Coverage

blib/lib/Math/PlanePath/HilbertCurve.pm

Criterion	Covered	Total	%
statement	102	119	85.7
branch	47	58	81.0
condition	6	10	60.0
subroutine	13	15	86.6
pod	5	5	100.0
total	173	207	83.5

# Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see .


# http://www.woollythoughts.com/afghans/peano.html
#     Knitting
#
# http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node36.html
#     Closed path, curved parts
#
# http://www.wolframalpha.com/entities/calculators/Peano_curve/jh/4o/im/
#     Curved corners tilted to a diamond, or is it an 8-step pattern?
#
# http://www.davidsalomon.name/DC2advertis/AppendC.pdf
#
# http://www.cut-the-knot.org/do_you_know/hilbert.shtml
#     Java applet
#

package Math::PlanePath::HilbertCurve;
use 5.004;
use strict;
#use List::Util 'max','min';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 129;
use Math::PlanePath;
use Math::PlanePath::Base::NSEW;
@ISA = ('Math::PlanePath::Base::NSEW',
        'Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow',
  'round_up_pow',
  'bit_split_lowtohigh',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';

# uncomment this to run the ### lines
#use Smart::Comments;


use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
*xy_is_visited = \&Math::PlanePath::Base::Generic::xy_is_visited_quad1;


#------------------------------------------------------------------------------

# state=0    3--2   plain
#               |
#            0--1
#
# state=4    1--2  transpose
#            |  |
#            0  3
#
# state=8
#
# state=12   3  0  rot180 + transpose
#            |  |
#            2--1
#
# generated by tools/hilbert-curve-table.pl
my @next_state = (4,0,0,12, 0,4,4,8, 12,8,8,4, 8,12,12,0);
my @digit_to_x = (0,1,1,0, 0,0,1,1, 1,0,0,1, 1,1,0,0);
my @digit_to_y = (0,0,1,1, 0,1,1,0, 1,1,0,0, 1,0,0,1);
my @yx_to_digit = (0,1,3,2, 0,3,1,2, 2,3,1,0, 2,1,3,0);
my @min_digit = (0,0,1,0, 0,1,3,2, 2,undef,undef,undef,
                 0,0,3,0, 0,2,1,1, 2,undef,undef,undef,
                 2,2,3,1, 0,0,1,0, 0,undef,undef,undef,
                 2,1,1,2, 0,0,3,0, 0);
my @max_digit = (0,1,1,3, 3,2,3,3, 2,undef,undef,undef,
                 0,3,3,1, 3,3,1,2, 2,undef,undef,undef,
                 2,3,3,2, 3,3,1,1, 0,undef,undef,undef,
                 2,2,1,3, 3,1,3,3, 0);

sub n_to_xy {
  my ($self, $n) = @_;
  ### HilbertCurve n_to_xy(): $n
  ### hex: sprintf "%#X", $n

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  my $int = int($n);
  $n -= $int;   # fraction part

  my @ndigits = digit_split_lowtohigh($int,4);
  my $state = ($#ndigits & 1 ? 4 : 0);
  my $dirstate   = ($#ndigits & 1 ? 0 : 4); # default if all $ndigit==3

  my (@xbits, @ybits);
  foreach my $i (reverse 0 .. $#ndigits) {    # digits high to low
    my $ndigit = $ndigits[$i];
    $state += $ndigit;
    if ($ndigit != 3) {
      $dirstate = $state;  # lowest non-3 digit
    }

    $xbits[$i] = $digit_to_x[$state];
    $ybits[$i] = $digit_to_y[$state];
    $state = $next_state[$state];
  }

  my $zero = ($int * 0); # inherit bigint 0
  return ($n * ($digit_to_x[$dirstate+1] - $digit_to_x[$dirstate]) # frac
          + digit_join_lowtohigh (\@xbits, 2, $zero),

          $n * ($digit_to_y[$dirstate+1] - $digit_to_y[$dirstate]) # frac
          + digit_join_lowtohigh (\@ybits, 2, $zero));
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### HilbertCurve xy_to_n(): "$x, $y"

  $x = round_nearest ($x);
  if (is_infinite($x)) { return $x; }
  $y = round_nearest ($y);
  if (is_infinite($y)) { return $y; }

  if ($x < 0 || $y < 0) {
    return undef;
  }

  my @xbits = bit_split_lowtohigh($x);
  my @ybits = bit_split_lowtohigh($y);
  my $numbits = max($#xbits,$#ybits);

  my @ndigits;
  my $state = ($numbits & 1 ? 4 : 0);
  foreach my $i (reverse 0 .. $numbits) {   # high to low
    ### at: "state=$state  xbit=".($xbits[$i]||0)." ybit=".($ybits[$i]||0)
    my $ndigit = $yx_to_digit[$state + 2*($ybits[$i]||0) + ($xbits[$i]||0)];
    $ndigits[$i] = $ndigit;
    $state = $next_state[$state+$ndigit];
  }
  ### @ndigits
  return digit_join_lowtohigh(\@ndigits, 4,
                              $x * 0 * $y); # inherit bignum 0
}


# rect_to_n_range() finds the exact minimum/maximum N in the given rectangle.
#
# The strategy is similar to xy_to_n(), except that at each bit position
# instead of taking a bit of x,y from the input instead those bits are
# chosen from among the 4 sub-parts according to which has the maximum N and
# is within the given target rectangle.  The final result is both an $n_max
# and a $x_max,$y_max which is its position, but only the $n_max is
# returned.
#
# At a given sub-part the comparisons ask whether x1 is above or below the
# midpoint, and likewise x2,y1,y2.  Since x2>=x1 and y2>=y1 there's only 3
# combinations of x1>=cmp,x2>=cmp, not 4.

# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### HilbertCurve rect_to_n_range(): "$x1,$y1, $x2,$y2"

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);
  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;

  if ($x2 < 0 || $y2 < 0) {
    return (1, 0); # rectangle outside first quadrant
  }

  my $n_min = my $n_max
    = my $x_min = my $y_min
      = my $x_max = my $y_max
        = ($x1 * 0 * $x2 * $y1 * $y2); # inherit bignum 0

  my ($len, $level) = round_down_pow (($x2 > $y2 ? $x2 : $y2),
                                      2);
  ### $len
  ### $level
  if (is_infinite($level)) {
    return (0, $level);
  }
  my $min_state = my $max_state = ($level & 1 ? 4 : 0);

  while ($level >= 0) {
    {
      my $x_cmp = $x_min + $len;
      my $y_cmp = $y_min + $len;
      my $digit = $min_digit[3*$min_state
                             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
                             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];

      $n_min = 4*$n_min + $digit;
      $min_state += $digit;
      if ($digit_to_x[$min_state]) { $x_min += $len; }
      if ($digit_to_y[$min_state]) { $y_min += $len; }
      $min_state = $next_state[$min_state];
    }
    {
      my $x_cmp = $x_max + $len;
      my $y_cmp = $y_max + $len;
      my $digit = $max_digit[3*$max_state
                             + ($x1 >= $x_cmp ? 2 : $x2 >= $x_cmp ? 1 : 0)
                             + ($y1 >= $y_cmp ? 6 : $y2 >= $y_cmp ? 3 : 0)];

      $n_max = 4*$n_max + $digit;
      $max_state += $digit;
      if ($digit_to_x[$max_state]) { $x_max += $len; }
      if ($digit_to_y[$max_state]) { $y_max += $len; }
      $max_state = $next_state[$max_state];
    }

    $len = int($len/2);
    $level--;
  }

  return ($n_min, $n_max);
}

#------------------------------------------------------------------------------

# shared by Math::PlanePath::AR2W2Curve and others
sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, 4**$level - 1);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  my ($pow, $exp) = round_up_pow ($n+1, 4);
  return $exp;
}

#------------------------------------------------------------------------------
1;
__END__