帰ってきたネコブログ

ジュリア集合を描いたフォルダが中間データのせいで2GBくらいあって、
ちょっとどうにかしたいのでスクリプトを書き直した。

use v5.14;
use strict;
use warnings;

use Imager;
use List::Util qw/max sum/;
use Time::HiRes qw/time/;

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;
    }

    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;

        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++) {
                NV z_r = (ix * d_r) + fP0_r;
                NV z_i = (iy * d_i) + fP0_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;
                }

                p[ix] = i;
            }

            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 = $_ / 255.0;
        Imager::Color->new( hue => 0.0, v => $val, s => 0.0 );
    } 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 => 2000;
use constant MARGIN => 10;
use constant DST_DIR => './eee';

# 出力先フォルダを生成
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 $n = 100;
my ( $w_mini, $h_mini ) = ( 250, 250 );
foreach ( 1..$n ) {
    my ( $a_r, $a_i ) = next_parameter();

    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}) < 32 ) {
        # 階調が少ないのでスキップ
        next;
    }

    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 );
    set_loop_limit( 1000 );

    my $src = generate( KS, KS, $a_r, $a_i );
    foreach my $curve ( 0.8, 1.0, 1.2 ) {
        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 );
    $i++;
}

完璧ではないにしろ、中間データはいらなくなった。
ログさえあれば、あとで何とでもなる。

というわけで、本日の1枚。

( re, im ) = ( 0.23, -0.532 )

2016/08/22 追記
生成したジュリア集合が画像サイズに合わせて良い感じに拡大されるように、
上記のスクリプトを修正しました。
— ここまで追記 —

おしまい。

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