piDay2022.groovy

#@ UIService ui
#@ DatasetService ds

/*
 * Happy Pi Day 2022
 * 
 * This visualizes the bifurcation diagram for the "standard circle map"
 * https://en.wikipedia.org/wiki/Arnold_tongue#Standard_circle_map
 * 
 * (and its parameters are expressed in terms of pi)
 * 
 * John Bogovic
 * 
 * This script was written for Fiji / ImageJ2.
 */
 
import net.imglib2.type.numeric.integer.*;
import net.imagej.axis.*;

import ij.*;
import ij.process.*;

import net.imglib2.img.array.*;
import net.imagej.*;

pi = Math.PI;
pi2 =  2 * pi;
nsamples = pi*pi*pi*pi*pi*pi as long;
nCircleSamples = pi*pi*pi*pi*pi*pi as long;
step = pi / (pi*pi*pi);
omega = pi / (pi*pi);
startK = pi / 2;
endK   = 4 * pi;

def circleMapOrbit( net.imglib2.RandomAccess ra, double omega, double K, long nsamples ) {
	
	t0 = 0.0;
	while( t0 < pi2 ) {
		t0 += step;
		tht = t0;
		(1..nsamples).each {
			tht = tht + omega + (( K / pi2 ) * Math.sin( pi2 * tht ));
			tht = tht % pi2;
			if (tht < 0 ) { 
				tht += pi2;
			}
			p = ((tht * nCircleSamples ) / pi2) as long;
			ra.setPosition( p, 0 );
			ra.get().inc();
		};
	}
}

img = ds.create(new ImgPlus( ArrayImgs.unsignedShorts( nCircleSamples, nsamples )));
ra = img.randomAccess();
for( int i = 0; i < nsamples; i++ ) {
	println ( "i: " + i );
	ra.setPosition( i, 1 );
	K = startK + ( endK - startK ) * i / nsamples;
	circleMapOrbit( ra, omega, K, nsamples );
}

ui.show( img );