Visual analogue of a Shepard tone

A Shepard tone is an auditory illusion that appears to indefinitely ascend or descend in pitch, without actually changing pitch at all.

 

Shepard tones work because they actually contain multiple tones, separated by octaves. As tones get higher in pitch, they fade out. New tones fade in at the lower pitches. The net effect is that it sounds like all the constituent tones are continually increasing in pitch -- and they are, but pitches fade in and out so that, on average, the pitch composition is constant.

Since 2D quasicrystals can be rendered as a sum of plane-waves, it is possible to form the analogue of a Shepard tone with these visual objects. Each plane wave is replaced with a collection of plane waves, at 2,4,8,16... etc times the spatial frequency of the original plane wave.

The relative amplitudes of the plane waves are set so that the spatial frequency stays approximately the same even as the underlying waves are scaled. The result is a quasicrystal that appears to zoom in or out indefinitely, without fundamentally changing in structure.

The infinite zoom effects creates a motion-fatigue optical illusion, which will cause illusory contraction of your visual field after staring at the GIF below:

 

 

 

More quasicrystal zoom GIFs can be found here. You can run and modify the code I used to generate these animation. Copy the following code into a file called QuasiZoom.java. Then, in a terminal, type "javac QuasiZoom.java" in the same directory, and then "java QuasiZoom". Various parameters to tune the output are noted in comments in the code. Then use Gimp to make an animated GIF.

import java.awt.Color;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import javax.imageio.ImageIO;
import static java.lang.Math.*;

public class QuasiZoom {

    // Defines a gaussian function. We will use this to define the
    // envelope of spatial frequencies
    public static double gaussian(double x) {
        return exp(-x*x/2)/sqrt(2*PI);
    }

    public static void main(String[] args) throws IOException {
        int k = 5;        //number of plane waves
        int stripes = 3;  //number of stripes per wave
        int N = 500;      //image size in pixels
        int divisions=40; //number of frames to divide the animation into
        int N2 = N/2;

        BufferedImage it = new BufferedImage(N, N, BufferedImage.TYPE_INT_RGB);

        //the range of different spatial frequencies
        int [] M=new int[]{1,2,4,8,16,32,64,128,256};
        
    //the main ( central ) spatial frequency
        double mean=log(16);

    //the spread of the spatial frequency envelope
        double sigma=1;

    //counts the frames 
        int ss=0;

    //iterate over spatial scales, scaling geometrically
        for (double sc=2.0; sc>1.0; sc/=pow(2,1./divisions)) 
        {    
            System.out.println("frame = "+ss);

            //adjust the  wavelengths for the current spatial scale
            double [] m=new double[M.length];
            for (int l=0; l<M.length; l++)
                m[l]=M[l]*sc;

            //modulate each wavelength by a gaussian envelop in log
            //frequency, centered around aforementioned mean with defined
            //standard deviation
            double sum=0;
            double [] W=new double[M.length];
            for (int l=0; l<M.length; l++) {
                W[l]=gaussian((log(m[l])-mean)/sigma);
                sum+=W[l];
            }
            sum*=k;

            for (int i = 0; i < N; i++) {
                for (int j = 0; j < N; j++) {

                    double x = j - N2, y = i - N2; //cartesian coordinates
                    double C = 0;                  // accumulator
 
                    // iterate over all k plane waves
                    for (double t = 0; t < PI; t += PI / k){
                        //compute the phase of the plane wave
                        double ph=(x*cos(t)+y*sin(t))*2*PI*stripes/N;
                        //take a weighted sum over the different spatial scales
                        for (int l=0; l<M.length; l++)
                            C += (cos(ph*m[l]))*W[l];
                    }
                    // convert the summed waves to a [0,1] interval
                    // and then convert to [0,255] greyscale color
                    C = min(1,max(0,(C*0.5+0.5)/sum));
                    int c = (int) (C * 255);
                    it.setRGB(i, j, c | (c << 8) | (c << 16));
                }
            }
            ImageIO.write(it, "png", new File("out"+(ss++)+".png"));
        }
        
    }
}

Animating quasicrystals

This post is inspired by the recent Nobel prize for the discovery of quasicrystals in nature.

We spoke briefly of quasicrystals previously : patterns that appear somewhat periodic but never actually repeat. Two-dimensional quasicrystals can be generated by summing together four or more plane waves. In collaboration with Keegan, we've explored what happens when you make these plane waves travel, and also viewing the quasicrystals in logarithmic coordinates as shown in this earlier post.

 

The Java code below generated sequences of images for animation. Once you have a series of images, you can use Gimp to make an animated GIF. (Convert to black-and-white and 8-color indexed mode before saving to reduce file size.)

import java.awt.image.BufferedImage;  
 import java.io.File;  
 import java.io.IOException;  
 import javax.imageio.ImageIO;  
 import static java.lang.Math.*;  
 public class QUASI1 {  
   public static void main(String [] args) throws IOException {  
     int k=4;     //numer of plane waves  
     int stripes = 27; //number of stripes per wave  
     int N = 800;   //image size in pixels  
     int N2 = N/2;  
     BufferedImage it = new BufferedImage(N,N,BufferedImage.TYPE_INT_RGB);  
     for (double phase=0; phase<2*PI; phase+=2*PI/30) {  
       for ( int i=0; i<N; i++ ) for ( int j=0; j<N; j++ ) {  
         double x = j-N2, y = i-N2; //cartesian coordinates  
         double theta = atan2(y,x); //log-polar coordinates  
         double r = log(sqrt(x*x+y*y));  
         double C=0;        // accumulator  
         for (double t=0; t<PI; t+=PI/k)  
           C+=cos((theta*cos(t)-r*sin(t))*stripes+phase);  
           // use the following line for cartesian crystals:  
           //C+=cos((x*cos(t)+y*sin(t))*2*PI*stripes/N+phase);  
         int c=(int)((C+k)/(k*2)*255);  
         it.setRGB(i,j,c|(c<<8)|(c<<16));  
       }  
       ImageIO.write(it,"png",new File("Test"+(int)(180*phase/PI)+".png")) ;  
     }  
   }  
 }  

To execute this code, copy it into a file names "QUASI1.java". Then, from the terminal, run "javac QUASI1.java", and finally "java QUASI1" to execute it. Or, paste it into your favorite Java IDE.

Changing K will change the degree of symmetry in the crystal. Changing N sets the size of the output images. Changing stripes sets how many wave cycles fit in the rendered image, with larger numbers leading to finer structure.

Since quasicrystals are aperiodic, it is not possible to wrap them around the log-polar "tunnel" ( seen on the left above ) such that two edges of the image meet perfectly. However, for quasicrystals composed of only a few plane waves, you can sometimes get two regions to align well enough to be unnoticeable, especially in black and white.

 

Be sure to check out Keegan's implementation, which doubles as a nice introduction to coding in Haskell.

Quasicrystal

Quasicrystals are aperiodic structures that still have some regularities. Both the spatial and frequency domain representations of quasicrystals are quite beautiful. See also this more recent post for more quasicrystal entertainment. Here is the classic diffraction pattern from a synthesized icosahedral quasicrystal

You can create two-dimensional quasicrystals by summing together more than 3 plane waves. For example, here are some spatial domain two-dimensional quasicrystals :

4-fold :
5-fold :
7-fold :

[overflow gallery]

The picture at the beginning of this post is an x-ray diffraction pattern, which is a little bit like the Fourier transform of a quasicrystal lattice. Here is a simulated 2D Fourier spectrum from a 7-fold 2D quasicrystal :

You can create your own images like this by rendering large two-dimension quasicrystals as described here or here. The frequency domain images can be rendered using ImageJ, or your own personal 2d FFT code. In ImageJ, open an image, and then create a frequency domain image by clicking Process→FFT→FFT. If you rendered a N-fold quasicrystal, you should see 2N points arranged in a circle around the origin, corresponding to your N plane waves. If you applied some nonlinear image operation, like contrast enhancement or thresholding, to the spatial domain image, you will have created some harmonics and overtones of your original N plane waves, which should appear as a constellation of other points that themselves are arranges on a quasi-crystal lattice.

Quasicrystals can be related to more abstract mathematical concepts. For instance, this article finds that some quasicrystals are related to sorting algoritms :

Current research in abstract tiling theory examines tilings of high rotational symmetry in collaboration with Remy Mosseri and co-workers. Possibly the limit of high rotational symmetry may prove easier to analyze than specific finite symmetries (10-fold, for example) of direct physical interest. Surprisingly, rhombus tilings are related to algorithms for sorting of lists. Counting the number of distinct tilings enumerates simultaneously the number of equivalance classes of sorting algorithms, a problem previously considered by computer scientist D.E. Knuth. Our random tiling theory implies an upper bound of log(2) for the tiling entropy per vertex, consistent with a conjecture by Knuth. Click here for a preprint on this research.