File Coverage

blib/lib/Math/DCT.pm

Criterion	Covered	Total	%
statement	52	52	100.0
branch	22	22	100.0
condition	16	18	88.8
subroutine	9	9	100.0
pod	5	5	100.0
total	104	106	98.1

package Math::DCT;

use 5.008;
use strict;
use warnings;

=head1 NAME

Math::DCT - 1D and NxN 2D Fast Discreet Cosine Transforms (DCT-II)

=head1 SYNOPSIS

    use Math::DCT qw/dct dct1d dct2d idct1d idct2d/;

    my $dct1d = dct([[1,2,3,4]]);
    $dct1d = dct1d([1,2,3,4]);

    my $dct2d = dct([[1,2],[3,4]]);
    $dct2d = dct2d([1,2,3,4]);

    my $idct1d = idct1d([1,2,3,4]);
    my $idct2d = idct2d([1,2,3,4]);

=head1 VERSION

Version 0.04_0

=cut

our $VERSION = '0.04_0';

require XSLoader;
XSLoader::load('Math::DCT', $VERSION);

use Exporter qw(import);

our @EXPORT_OK = qw(
    dct
    dct1d
    dct2d
    idct1d
    idct2d
);

our %EXPORT_TAGS;
$EXPORT_TAGS{all} = [@EXPORT_OK];

=head1 DESCRIPTION

An unscaled DCT-II implementation for 1D and NxN 2D matrices implemented in XS.
For array sizes which are a power of 2, a fast algorithm (FCT) described by Lee is
used (with the addition of a coefficient table that makes it even faster than some
common implementations), plus an unscaled version of the Arai, Agui, Nakajima
algorithm is used for size 1x8, 8x8 matrices. A less optimized algorithm is used
for the generic case, so any 1D or square 2D matrix can be processed.

For convenience the inverse functions are provided (inverse DCT-II, usually called
iDCT, is essentially a scaled DCT-III), although using just one generic algorithm.

The module was written for a perceptual hash project that needed 32x32 DCT-II, and
on a 2.5GHz i7 2015 Macbook Pro about 18000/s per thread are processed (C function).
The common 8x8 DCT-II uses a special path, (about 380000/s on that same CPU),
although for most image/video applications that require 8x8 DCT there are much faster
implementations (SIMD, approximations etc) that usually produce an already scaled result
for the specific application.

None of the algorithms used on this module are approximate, the test suite verifies
against a naive DCT-II implementation with a tolerance of 1e-08.

=head1 METHODS

=head2 C

  my $dct = dct([[1,2],[3,4]]);

Pass an array (ref) of either a single array, or N N-length arrays for 1D and 2D
DCT-II calculation respectivelly. The output will be an arrayref of array(s) with
the result of the transform.

=cut

sub dct {
    my $vector = shift;
    die "Expect array of array(s)"
        unless $vector && ref $vector eq 'ARRAY'
        && $vector->[0] && ref $vector->[0] eq 'ARRAY';

    my $dim = scalar(@$vector);
    my $sz  = scalar(@{$vector->[0]});
    die "Expect 1d or NxN 2d arrays" unless $dim == 1 || $dim == $sz;

    my $pack;
    for (my $x = 0; $x < $dim; $x++) {
        $pack .= pack "d$sz", @{$vector->[$x]}
    }

    if ($dim > 1) {
        $sz == 8
            ? fct8_2d($pack)
            : 0 == ($sz & ($sz - 1)) ? fast_dct_2d($pack, $sz) : dct_2d($pack, $sz);
    } else {
        $sz == 8
            ? fct8_1d($pack)
            : 0 == ($sz & ($sz - 1)) ? fast_dct_1d($pack, $sz) : dct_1d($pack, $sz);
    }

    my $result;
    foreach (0..$dim-1) {
        $result->[$_] = [unpack "d".($sz), substr $pack, $_ * $sz*8, $sz*8];
    }
    return $result;
}

=head2 C

  my $dct = dct1d([1,2,3]);

Pass an array (ref) for a 1D DCT-II calculation. The output will be an arrayref
with the result of the transform.

=cut

sub dct1d {
    my $input = shift;
    my $sz    = scalar @$input;
    my $pack  = pack "d$sz", @$input;
    $sz == 8
        ? fct8_1d($pack)
        : 0 == ($sz & ($sz - 1)) ? fast_dct_1d($pack, $sz) : dct_1d($pack, $sz);
    my @result = unpack "d$sz", $pack;
    return \@result;
}

=head2 C

  my $idct = idct1d([1,2,3]);

Pass an array (ref) for a 1D iDCT-II calculation. The output will be an arrayref
with the result of the transform. This is essentially a DCT-III transform scaled
by 2/N.

=cut

sub idct1d {
    my $input = shift;
    my $sz    = scalar @$input;
    my $pack  = pack "d$sz", @$input;
    idct_1d($pack, $sz);
    my @result = unpack "d$sz", $pack;
    return \@result;
}

=head2 C

  my $dct = dct2d(
      [1,2,3,4],   # Arrayref containing your NxN matrix
      2            # Optionally, the size N of your array (sqrt of its length)
  );

Pass an array (ref) for a 2D DCT-II calculation. The length of the array is expected
to be a square (as only NxN arrays are supported) - you can optionally pass N as
the second argument to avoid a C calculation.
The output will be an arrayref with the result of the transform.

If your 2D data is available in a 1D array as is usual with most image manipulation
etc cases, this function will be faster than C, as the DCT calculation is
in any case done on a 1D array, hence you save the cost of conversion.

=cut

sub dct2d {
    my $input = shift;
    my $sz    = shift || sqrt(scalar @$input);
    my $pack  = pack "d".($sz*$sz), @$input;
    $sz == 8
        ? fct8_2d($pack)
        : 0 == ($sz & ($sz - 1)) ? fast_dct_2d($pack, $sz) : dct_2d($pack, $sz);
    my @result = unpack "d".($sz*$sz), $pack;
    return \@result;
}

=head2 C

  my $idct = idct2d(
      [1,2,3,4],   # Arrayref containing your NxN matrix
      2            # Optionally, the size N of your array (sqrt of its length)
  );

Pass an array (ref) for a 2D iDCT-II calculation. The length of the array is expected
to be a square (as only NxN arrays are supported) - you can optionally pass N as
the second argument to avoid a C calculation.
This is essentially a DCT-III transform scaled by 2/N.
The output will be an arrayref with the result of the transform.

=cut

sub idct2d {
    my $input = shift;
    my $sz    = shift || sqrt(scalar @$input);
    my $pack  = pack "d".($sz*$sz), @$input;
    idct_2d($pack, $sz);
    my @result = unpack "d".($sz*$sz), $pack;
    return \@result;
}

=head1 USAGE NOTES

The C functions are not exported, but theoretically you could use them directly
if you do your own C. The fast versions for power-of-2 size arrays
are C and C, while the generic versions are C
and C (with inverse functions as C and C).
The specialized size-8 versions are C and C.
First argument is a C (use C), second is the size N (except
for the fct8* functions which don't need a second argument).

=head1 ACKNOWLEDGEMENTS

C code for 1D DCT from Project Nayuki was adapted and improved where possible.

(L)

=head1 AUTHOR

Dimitrios Kechagias, C<<  >>

=head1 BUGS

Please report any bugs or feature requests to C,
or through the web interface at L.
I will be notified, and then you'll automatically be notified of progress on your
bug as I make changes.

=head1 GIT

L

=head1 COPYRIGHT & LICENSE

Copyright (C) 2019, SpareRoom.com

This program is free software; you can redistribute
it and/or modify it under the same terms as Perl itself.

=cut

1;