帰ってきたネコブログ

ジュリア集合もアンチエイリアシングして滑らかにしようと思ったんだけど、
そもそもギザギザの部分が延々に繰り返すわけで、
結果的にすべてギザギザになるのは正しいってことに気づいて、
とりあえず着色してみたら気にならなくなりました。

前回は、大きな画像に描画して、縮小処理時に滑らかにする作戦でしたが、
今回は周辺の計算を行って平均をとることにしました。

use v5.14;
use strict;
use warnings;

use Imager;
use List::Util qw/max sum/;
use Time::HiRes qw/time/;
use Getopt::Long qw/:config posix_default no_ignore_case bundling auto_help/;

use Inline C => Config => LIBS => '-lm' => OPTIMIZE => '-O';
use Inline C => q{
    #include <stdio.h>
    #include <math.h>

    static IV fLoopLimit = 1000;
    static NV fP0_r = 0.0; // 実部（x方向）
    static NV fP0_i = 0.0; // 虚部（y方向）
    static NV fP1_r = 0.5;
    static NV fP1_i = 0.5;

    void set_loop_limit(IV limit)
    {
        fLoopLimit = limit;
    }

    void set_render_area(NV p0_r, NV p0_i, NV p1_r, NV p1_i)
    {
        fP0_r = p0_r;
        fP0_i = p0_i;
        fP1_r = p1_r;
        fP1_i = p1_i;
    }

    IV calculate(NV z_r, NV z_i, NV a_r, NV a_i)
    {
        IV i = 0;
        for (; i<fLoopLimit; i++) {
            const NV z2_r = (z_r * z_r) - (z_i * z_i) + a_r;
            const NV z2_i = (2.0 * z_r * z_i) + a_i;

            if ( 4.0 < ((z2_r * z2_r) + (z2_i * z2_i)) ) {
                break;
            }

            z_r = z2_r, z_i = z2_i;
        }

        return i;
    }

    SV * generate(IV w, IV h, NV a_r, NV a_i)
    {
        const NV d_r = ( fP1_r - fP0_r ) / w;
        const NV d_i = ( fP1_i - fP0_i ) / h;
        const NV dd_r = d_r / 3.0;
        const NV dd_i = d_i / 3.0;
        const NV dd2_r = dd_r / sqrt(2.0);
        const NV dd2_i = dd_i / sqrt(2.0);

        IV *dst = NULL;
        Newx( dst, (w * h), IV );

        printf( "progress: %4d/%4d", 0, h );
        fflush( stdout );

        for (IV iy=0; iy<h; iy++) {
            IV *p = dst + (w * iy);
            for (IV ix=0; ix<h; ix++) {
                const NV z_r = (ix * d_r) + fP0_r;
                const NV z_i = (iy * d_i) + fP0_i;

                IV wk = 0;
                wk = calculate( z_r, z_i, a_r, a_i ) * 8.0;
                wk += calculate( z_r + dd_r, z_i, a_r, a_i );
                wk += calculate( z_r - dd_r, z_i, a_r, a_i );
                wk += calculate( z_r, z_i + dd_i, a_r, a_i );
                wk += calculate( z_r, z_i - dd_i, a_r, a_i );                
                wk += calculate( z_r + dd2_r, z_i + dd2_i, a_r, a_i );
                wk += calculate( z_r + dd2_r, z_i - dd2_i, a_r, a_i );
                wk += calculate( z_r - dd2_r, z_i + dd2_i, a_r, a_i );
                wk += calculate( z_r - dd2_r, z_i - dd2_i, a_r, a_i );

                p[ix] = wk;
            }

            printf( "\rprogress: %4d/%4d", iy + 1, h );
            fflush( stdout );
        }

        printf( " -> calclated!\n" );

        SV *ret = newAV();
        av_fill( ret, (w * h) - 1 );
        SV **base = AvARRAY( ret );
        for (IV iy=0; iy<h; iy++) {
            const IV *p = dst + (w * iy);
            for (IV ix=0; ix<w; ix++) {
                (*base++) = newSViv( p[ix] );
            }
        }

        Safefree( dst );
        return newRV_noinc( (SV *)ret );
    }
};

sub render_image {
    my ( $src, $w, $h, $th, $curve, $palette ) = @_;

    my $img = Imager->new(
        xsize => $w, ysize => $h, channels => 3 );

    my $n = scalar( @{$palette} );
    my @dst = map {
        my $i = int( (($_ / $th) ** $curve) * ($n - 1) );
        $palette->[ ($n <= $i) ? ($n - 1) : $i ];
    } @{$src};

    my $iy = 0;
    while ( my @line = splice(@dst, 0, $w) ) {
        $img->setscanline( y => $iy, pixels => \@line );
        $iy++;
    }

    return $img;
}

sub save_image {
    my ( $dst_file, $src, $w, $h, $hist, $curve ) = @_;

    my $th = scalar(@{$hist}) - 1;
    $th-- while ( not $hist->[$th] );

    my @palette = map {
        my $val = ($_ < 150) ? ($_ / 150.0) : 1.0;
        my $hue = -6.0 + (($_ / 255.0) * 80);
        my $sat = ($_ < 180) ? 0.5 : 0.5 - (($_ - 180) / 255.0);
        Imager::Color->new( hue => $hue, v => $val, s => $sat );
    } 0..255;

    my $img = render_image( $src, $w, $h, $th, $curve, \@palette );
    $img->write( file => $dst_file ) or die $img->errstr;
    say 'wrote: ', $dst_file;
}

sub calc_fit_size {
    my ( $src, $w, $h ) = @_;

    # 端っこは、背景色と仮定
    my $ignore = $src->[$w - 1];

    my ( $x_min, $y_min, $x_max, $y_max ) = ( $w / 2, $h / 2, $w / 2, $h / 2 );
    for (my $iy=0; $iy<$h; $iy++) {
        my $offset = $iy * $w;

        my $found = 0;
        for (my $ix=0; $ix<$w; $ix++) {
            my $val = $src->[$offset + $ix];
            next if $val == $ignore;

            $x_min = $ix if $ix < $x_min;
            $x_max = $ix if $x_max < $ix;
            $found = 1;
        }

        if ( $found ) {
            $y_min = $iy if $iy < $y_min;
            $y_max = $iy if $y_max < $iy;
        }
    }

    return ( $x_max - $x_min + 1, $y_max - $y_min + 1 );
}

sub get_hist {
    my ( $src ) = @_;

    my $val_max = int( max(@{$src}) );
    my @hist = map { 0; } 0..$val_max;
    $hist[int($_)]++ for @{$src};

    return \@hist;
}

sub append_log {
    my ( $log_file, $log ) = @_;

    open( my $fh, '>>', $log_file ) or die "cannot open $log_file : $!";
    print $fh $log, "\n";
    close( $fh );
}

sub get_latest_log_no {
    my $log_file = shift;

    return -1 if ( not -e $log_file );

    local $/ = "\n";
    open( my $fh, '<', $log_file ) or die "cannot open $log_file : $!";
    my @logs = <$fh>;
    close( $fh );

    chomp @logs;
    @logs = grep { $_; } @logs;
    if ( not @logs ) {
        return -1;
    }
    else {
        my @tmp = split( /,/, pop(@logs) );
        return int( shift(@tmp) );
    }
}

sub next_parameter {
    # 生成ログが見やすいように任意の精度で切り捨て
    return map {
        sprintf('%.3f', (2 * rand() - 1.0)) * 1.0;
    } 1..2;
}

use constant KS => 640;
use constant MARGIN => 10;

GetOptions(
    \my %opt,
    'area=s',
    'param=s',
    'dst=s'
);

if ( not exists $opt{dst} ) {
    die '"dst" is required.';
}

my ( $a_r, $a_i ) = ( exists $opt{param} )
    ? split(/,/, $opt{param}) : next_parameter();
my $dst_dir = $opt{dst};

say 'param = ', join(',', $a_r, $a_i);

# 出力先フォルダを生成
if ( not -e $dst_dir ) {
    mkdir $dst_dir;
}

local $| = 1;
my $log_file = join( '/', $dst_dir, 'log.txt' );
my $i = get_latest_log_no($log_file) + 1;
my ( $w_mini, $h_mini ) = ( 250, 250 );

my $r = 2.0;
set_render_area( -$r, -$r, $r, $r );
set_loop_limit( 500 );

my $src_mini = generate( $w_mini, $h_mini, $a_r, $a_i );
my $hist = get_hist( $src_mini );
if ( (grep { 0 < $_; } @{$hist}) < 2 ) {
    die "cannot render!";
}

if ( not exists $opt{area} ) {
    my ( $w_fit, $h_fit ) = calc_fit_size( $src_mini, $w_mini, $h_mini );
    $r = max( $r * ($w_fit / $w_mini), $r * ($h_fit / $h_mini) );

    # ぎりぎりではなく、余白ができるように調整
    $r = abs( ((KS - MARGIN) / KS) * $r );
    $r = sprintf( '%.3f', $r );

    say "area : -$r, -$r, $r, $r";
    set_render_area( -$r, -$r, $r, $r );
}
else {
    my @tmp = split( /,/, $opt{area} );
    die '"area" is invalid.' if scalar(@tmp) != 4;

    set_render_area( @tmp );
}

set_loop_limit( 5000 );
my $src = generate( KS, KS, $a_r, $a_i );
foreach my $curve ( 1.2, 1.4, 1.6 ) {
    my $dst_file = sprintf( '%03d_%.1f.png', $i, $curve );
    $dst_file = join( '/', $dst_dir, $dst_file );
    save_image( $dst_file, $src, KS, KS, $hist, $curve );
}

my $log = join( ',', sprintf('%03d', $i), $a_r, $a_i );
append_log( $log_file, $log );

一応、描画範囲の指定もできますがGUIじゃないと面倒ですね。
あと、周辺の計算を行って平均をとった結果、
大きな画像を作って縮小処理をした感じになったので、
ギザギザが少し緩和された気がします。

あと、着色の方はありがちな感じごめんなさい。
今回から出力先ディレクトリの指定が必須です。
$ perl aaa.pl --param="-0.034,-0.657" --dst="test"

おしまい。

カテゴリー: Perl, コンピューター by ネコ on 2016年08月22日 | タグ: Imager > フラクタル > メモ