Representation of Maximally Regular Textures in Human Visual Cortex

Wallpaper stimulus generation.

Wallpaper patterns are repetitive 2D patterns that tile the plane. There are 17 unique wallpaper patterns, corresponding to the set of planar symmetry groups (Fedorov, 1891; Polya, 1924; Liu et al., 2010). Each of the 17 wallpaper groups is built on one of five types of “unit lattices” that are used for tiling the plane without gaps. The “fundamental region” is the smallest repeating region in the wallpaper patterns. Within the unit lattice, multiple rigid transformations are applied to the fundamental region, which give rise to symmetries within the wallpaper group. As mentioned earlier, each wallpaper group contains a distinct combination of the following four fundamental symmetries: translations, rotations, reflections, and glide reflections.

In the experiments presented here, we measured neural responses to a subset of four of the 17 wallpaper groups: P2, P3, P4, and P6. These four groups all contain rotation symmetries, as well as translation symmetries given by the tiling of the unit lattice, but no other symmetries. Rotation symmetry around a particular point can be defined in terms of its order n, which means that the fundamental region can be rotated by an angle 360°/n without changing. The four wallpaper groups used in the present experiments are illustrated schematically in Figure 1, using a comma-shaped symbol as the fundamental region. Each group contains rotation symmetry around several points that vary in order. For P2, the maximum order of rotation symmetry is 2; for P3, it is 3; for P4, it is 4; and for P6, it is 6.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Tutorial examples of the four distinct wallpaper patterns used in the experiment, containing no reflection or glide symmetries. The maximum order of rotation symmetry for each wallpaper group is indicated next to each image. In these examples, the fundamental region is a comma-like symbol (illustration adapted from Wade, 1993).

We generated exemplars from the four wallpaper groups by using random noise textures as the fundamental regions based on a modification of the stimuli developed by Clarke et al. (2011b). In Figure 2, we show an exemplar from each of the four wallpaper groups, and indicate the fundamental region, rotation symmetry centers, and unit lattice on a magnified region of each texture. Although all four exemplars are generated from a fundamental region containing random noise, the resulting wallpapers are quite visually appealing, and very different from one another. Importantly, using random noise in the fundamental region means that for each wallpaper group an almost infinite number of distinct exemplars can be easily generated (a range of exemplars from group P6 are shown in Fig. 4).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

The stimuli used in the experiments reported here. Bottom, An exemplar from each of the four wallpaper groups. Top, Part of each exemplar, magnified so that the details of the pattern are easier to see. The fundamental regions are indicated with colored shading. The rotation symmetries in each pattern are indicated as follows: rotation of order 2 (rhombus), rotation of order 3 (triangle), rotation of order 4 (square), and rotation of order 6 (hexagon). The unit lattice is indicated with a yellow outline. Note how the fundamental region is rotated and repeated within the unit lattice, which is then used to tile the plane in each of the four exemplars. In the case of P3, the fundamental region is cut in half, and each half is rotated and repeated separately.

Generating wallpaper patterns by directly tiling the unit lattice is algorithmically challenging because unit lattices can have non-right angles, making them difficult to concatenate on a raster graphics display. We circumvented this issue by generating the smallest rectangular tile that contains the unit lattice and that can be easily replicated, here referred to as the “repeating tile.” The generation of a repeating tile starts with the fundamental region, which is cut from a larger “base texture.” The base texture is generated by convolving spatial white noise with a Gaussian-like filter constrained by the dimensions (width and height) of the base texture.

The fundamental region forms the basis for the unit lattice, which in turn constrains the repeating tile. Because of this, the shape and size of the fundamental region depends on the wallpaper group being portrayed as well as the area of the repeating tile. The generation of the fundamental region for group P6, an isosceles triangle with angles 30°, 30°, and 120°, is shown in Step 1 of Figure 3. After the fundamental region is generated, it is then rotated and repeated to generate the repeating tile appropriate for a specific wallpaper group. This process is illustrated in Steps 2–5 of Figure 3. The finished repeating tile is then used to tile the plane (Fig. 3, Step 6), before a final postprocessing procedure is applied (Fig. 3, Step 7), as described below.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Diagram of stimulus generation steps for wallpaper group P6, as follows: 1, the fundamental region is selected from the filtered random-noise base texture; 2a, the fundamental region is rotated 120° and 240°, and combined with itself into an equilateral triangle; 2b, building the triangle “liner.” The fundamental region is rotated 180°; the top half is stacked at the bottom, and the bottom half goes to the top. This set is now combined with 300° and 60° rotated versions of the fundamental region. 3, Triangle and triangle liner are joined together to produce half of the repeating tile; 4, half of the repeating tile is rotated 180° and combined with itself to produce a rectangular repeating tile; 5, the finished repeating tile; 6, the repeating tile is used to tile the plane, repeated eight times down and five times across; and 7, the finished stimulus after postprocessing (power spectrum replacement, low-pass filtering, histogram equalization, and application of a circular aperture mask).

The dimensions of the fundamental region define the area of the repeating tile, which in turn defines the spatial scale of the exemplar. A key issue when designing the wallpapers was to parameterize the dimensions of the fundamental region such that spatial scale was uniform across different groups. Because of forced relationships among the area of the unit lattice, the area of the fundamental region, and the dimensions of the repeating tile, it is not possible to equate all parameters across all groups. For our stimuli, we decided to equalize the area of the repeating tile across groups to 100 × 100 image pixels. This meant that the unit lattices had an area of 100² pixels for the P2 and P4 groups, and 2/3 × 100² for the P3 and P6 groups (Fig. 2). The initial horizontal and vertical dimensions of the base texture were set to the square root of the tile area for all groups to ensure that the granularity within the textures was consistent. We calculated the area of the repeating tile as a function of the dimensions of the fundamental region for each group and cropped the base texture so that the repeating tile would have the desired area. For the wallpaper groups containing rotations of order 3 and 6 (P3 and P6 groups), an additional upscaling operation was performed: after cropping, the base texture was scaled up multiple times (usually 9–10 times), allowing us to concatenate the results of rotation with high precision (within 3 pixels). The finished repeating tile was then scaled back down so that the area of the repeating tile was matched across all groups.

For each of the four wallpaper groups, we generated 20 exemplars, each starting with a new randomly generated base texture and resulting fundamental region. We then applied the following postprocessing procedure to each of the exemplars to accomplish several goals. First, we performed a Fourier decomposition on each of the 20 exemplar images and computed the average power spectrum across all exemplars within each group. We then replaced the power spectrum of each individual exemplar with the average, so that the power spectrum was equated across all exemplars. To minimize edge artifacts due to tiling, each exemplar image was then low-pass filtered using a Gaussian filter with size 9 × 9 pixels and SD σ = 1 pixel. Using the histeq function in MATLAB, we performed histogram equalization to increase the visual salience of the patterns, and then scaled the gray-level pixel intensities to 2–255. Finally, we applied a circular aperture mask to the images to minimize interactions between the wallpaper group lattices and the pattern edges. The effect of the postprocessing procedure can be observed by comparing Steps 6 and 7 in Figure 3.

Control stimulus.

Phase-randomized control exemplars were generated for each of the 20 exemplars from each wallpaper group. Control exemplars had the same power spectrum as the exemplar images for each group. Although the phase scrambling eliminates the symmetry relationships within the unit lattice, the periodicity of the lattice itself is inherited by the phase-scrambled images. This means that all control exemplars, regardless of which wallpaper group they are derived from, degenerate to another symmetry group, namely P1. P1 is the simplest of the wallpaper groups and contains only translations of a region whose shape derives from the lattice. Because the different wallpaper groups have different lattices, these P1 controls are slightly different between groups. Among the four groups we studied here, P2 and P4 shared one lattice, and P3 and P6 shared another, different lattice. As we shall see, our experimental design takes these lattice differences into account by comparing the neural responses evoked by a wallpaper group to that evoked by the matched control exemplars that share the same lattice structure, but have no rotation symmetry.

Stimulus presentation.

In the functional MRI experiment, the wallpaper images were shown on a 47 inch Resonance Technology LCD display. Participants viewed the screen through a mirror at a distance of 277 cm. The screen had a resolution of 1024 × 768 pixels, an 8 bit color depth, and a refresh rate of 60 Hz. The mean luminance was 34.39 cd/m², and contrast was 90%. The circular aperture defining the edge of the wallpapers had a diameter of 11.2° of visual angle, and the rest of the screen was uniformly gray. In the EEG experiment, the stimuli were shown on a 24.5 inch Sony Trimaster EL PVM-2541 organic light-emitting diode display, with a screen resolution of 1920 × 1080 pixels, 8-bit color depth and a refresh rate of 60 Hz, viewed at a distance of 70 cm. The mean luminance was 69.93 cd/m², and contrast was 95%. The diameter of the circular aperture was slightly larger in this experiment at 13.8° of visual angle, and again the background was gray.

Structural and functional MRI acquisition.

Functional and structural MRI data were collected on a Discovery 750 MRI system (General Electric Healthcare) equipped with a 32-channel head coil (Nova Medical) at the Center for Cognitive and Neurobiological Imaging at Stanford University. For each participant, we acquired two whole-brain T1-weighted structural datasets (1.0 × 1.0 × 1.0 mm resolution, TE = 2.5 ms, TR = 6.6 ms, flip angle = 12°, FOV = 256 × 256 mm) that were used for tissue segmentation and registration with the functional scans. We also acquired a single whole-brain, T2-weighted structural dataset (1.0 × 1.0 × 1.0 mm resolution, TE = 75 ms, TR = 2500 ms, flip angle = 90°, FOV = 256 × 256 mm), which was used for defining the tissue boundaries used for the EEG head model. For the functional MRI experiment, we collected functional data consisting of 30 coronal slices (2.0 × 2.0 × 2.0 mm resolution, TE = 30 ms, TR = 2000 ms, flip angle = 77°, FOV = 220 × 220 mm), positioned at an oblique angle to maximize the coverage of occipital, ventral, and parietal cortices. The scans used for localization of functionally defined regions of interest (ROIs) had the same voxel size and TE/TR as that used in the main experiment, but a multiplexed EPI sequence was used that allowed for whole-brain coverage.

fMRI experimental procedure.

We used a block design for the fMRI experiment in which 12 s experiment blocks (referred to as “A blocks”) alternated periodically with 12 s control blocks (referred to as “B blocks”), yielding a 24 s base period for the paradigm, which was repeated 10 times in what we refer to as a “scan.” The paradigm is illustrated schematically in Figure 4. Ten stimulus cycles were shown per scan, with an additional half-cycle (one 12 s control block) being shown in the beginning of the scan to allow the brain and the scanner to settle. The data collected during this “dummy” period were not included in the analysis. Fourier analysis was used to quantify BOLD activation at frequencies equal to exact integer multiples of the 240 s scan period (0.0042 Hz) up to 0.5 Hz, as limited by the 2 s TR. The paradigm frequency was 10 cycles per scan, and all plots use this convention for the frequency axis.

During the experiment blocks, 10 wallpaper group exemplars were shown interleaved with their matched control exemplars, alternating at 1 Hz (e.g., 500 ms duration P1 control images followed by 500 ms presentation of the corresponding intact exemplar). During the control blocks, P1 control exemplars were alternated with another set of P1 control exemplars generated using the same procedure (for a total of 20 control images per block, with each image again presented for 500 ms). Experiment and control blocks differed only in the presence of intact rotation symmetry in the experiment block. Ten image pairs (control–exemplar or control–control) were shown in each block, always in the same order, and the first and last pairs were shown twice, for a total of 24 image presentations within a block. Within each scan, exemplars and control stimuli all came from the same wallpaper group, and we collected three 4 min scans for each group. To avoid the effects of scan order within a session, we scanned all four groups in a specific order, which was repeated three times. We always used the same repeating order (P2 → P3 → P4 → P6), but different starting groups were used for each participant, such that all four groups were used as the starting group for equal amounts of time. The fMRI experiment design is illustrated in Figure 4B.

We used a concurrent task to ensure that participants were alert and paying attention to the visual stimulus. Two times per block, for both A and B blocks, a randomly chosen image pair was shown at reduced contrast, and participants were instructed to press a button whenever they detected a contrast reduction. Before each scan, we adjusted the contrast reduction such that the average value for participant accuracy was kept at ∼85% correct.

EEG experimental procedure.

The EEG experiment was designed to mimic the fMRI experiment as closely as possible, while taking into account the superior temporal resolution of EEG. We used a steady-state design, in which P1 control images alternated with test images from each of the four wallpaper groups under study (P2, P3, P4, and P6). As in the fMRI study, each exemplar image was always preceded by its matched P1 control image, which was generated as described above. Each image was shown for 600 ms, and a control image followed by an exemplar image constituted a single stimulus cycle, with a frequency of 0.83 Hz. A trial consisted of 10 stimulus cycles corresponding to 10 different exemplar images and matched controls, which all came from the same wallpaper group. The images were always shown in the same order within a trial.

Participants initiated each trial with a button press, which allowed them to take breaks between trials. Trials from a single wallpaper group were presented in blocks of four repetitions, which were themselves repeated twice per session, and were shown in random order within each session. During each trial, participants performed the same concurrent task as in the fMRI experiment, as follows: two times per trial, an image pair was shown at reduced contrast, and the participants were instructed to press a button on the joypad. As in the fMRI experiment, we adjusted the contrast reduction such that the average value for participant accuracy was kept at ∼85% correct.

fMRI analysis.

After removing the dummy TR values, the fMRI data were preprocessed in AFNI (Cox, 1996). Preprocessing included the following steps: slice time correction; motion registration (the third TR of the first scan was always used as the base); scaling (each voxel was scaled to a mean of 100, and values were clipped at 200); and detrending [removing components corresponding to the six motion registration parameters, as well as Legendre polynomials of order 0 (constant signal), 1 (linear drift), and 2].

ROIs were functionally defined for each participant based on data from separate scanning sessions and registered to the experimental data by aligning the ROI anatomy with the experimental data using a rigid body transformation, which was then applied to the ROIs so that they were in registration with the functional data collected in the symmetry fMRI experiment.

The remainder of the analysis was performed in MATLAB. The time course data were first averaged across the three scans for each condition, and then across the voxels within each ROI. We then applied the Fourier transform to the average time course for each ROI, omitted DC, multiplied the amplitudes by 2 to get the single-sided spectrum, and scaled by dividing the amplitudes with the number of samples in the time course. An example of the resulting amplitude spectrum is shown in Figure 5. After computing the amplitude spectrum for all four conditions, we took the amplitude of two sidebands on both sides of the stimulus frequency over all four conditions, and averaging these 16 values (four sidebands × four conditions) to get a single number representing the noise level across all conditions. We then computed the signal-to-noise ratio (SNR) for each condition by dividing the amplitude at the stimulus frequency with the noise level. We performed this procedure for each of the 12 participants to generate participant-specific SNR values for all four conditions.

Functionally defined visual regions of interest.

Nine visual ROIs were defined individually for each of the 12 participants who took part in the fMRI experiments. The majority of the areas were defined based on retinotopic mapping, which was performed using the population receptive field method (Dumoulin and Wandell, 2008), with the exception of two of our participants for whom the traditional phase-mapping method with rings and wedges was used (DeYoe et al., 1994; Engel et al., 1994). The areas were defined manually based on standard criteria (Sereno et al., 1995; Engel et al., 1997): contralateral quarterfield representation for V1d, V1v, V2d, V2v, V3d, and V3v; contralateral hemifield representation for V4 (Wade et al., 2002); and V3A/B (Tootell et al., 1997). IPS0 was defined as an additional contralateral hemifield representation anterior to V3A/B (IPS0 was originally known as V7; Tootell et al., 1998; Wandell et al., 2007). VO1 as an additional contralateral hemifield representation abutting V4 (Brewer et al., 2005). We note that in three of our participants, IPS0 and/or VO1 could not be defined bilaterally, because of variability in the quality of the acquired data. Every ROI could be defined bilaterally in at least 10 participants.

In addition to the retinotopic areas, we defined LOC and motion-sensitive middle temporal cortex (hMT⁺) using fMRI block designs. During the LOC localizer scans, participants viewed blocks of images depicting common objects (12 s/block), alternating with blocks containing scrambled versions of the same objects. The stimuli were identical to the ones used in a previous study (Kourtzi and Kanwisher, 2000). hMT⁺ was defined using stimuli similar to those described by Huk and colleagues (2002): localizer scans blocks of low-contrast moving dots alternated with blocks of low-contrast static dots. Three repetitions of each of the two types of localizer scans were performed for each participant, which was sufficient to define LOC and hMT⁺ bilaterally in all participants.

EEG acquisition and preprocessing.

The EEG data were collected with 128-sensor HydroCell Sensor Nets (Electrical Geodesics) and bandpass filtered from 0.3 to 50 Hz. Following each experimental session, the 3D locations of all electrodes and three major fiducials (nasion, and left and right preauricular points) were digitized using a 3Space Fastrack 3D digitizer (Polhemus). The digitized locations were used to construct the EEG forward model (see EEG source-imaging analysis section). Raw data were subjected to an off-line sample-by-sample thresholding procedure in which noisy sensors were replaced by the average of the six nearest spatial neighbors. On average, <5% of the electrodes were substituted; these electrodes were mainly located near the forehead or the ears, and, as such, are likely to have a negligible impact on our results, as our stimuli will likely drive responses mainly at electrodes over occipital, temporal, and parietal locations. The EEG data were then re-referenced to the common average of all the sensors and segmented into 10 1.2-s-long epochs (each corresponding to exactly one cycle of image modulation). Epochs for which a large percentage of data samples exceeded a noise threshold (depending on the participant and ranging between 25 and 50 μV) were excluded from the analysis on a sensor-by-sensor basis. This was typically the case for epochs containing artifacts, such as blinks or eye movements.

A Fourier analysis was applied on every epoch using a discrete Fourier transform with a rectangular window. Given that each epoch was 1.2 s, this Fourier transformation led to a frequency resolution of δf = 0.42 Hz. For each frequency bin, the complex-valued Fourier coefficients were then averaged across all the epochs and all the trials. Thus, these average Fourier coefficients for each session were obtained from up to 80 data samples (10 epochs × 8 trials). Each participant did two sessions, and the Fourier coefficients from each session were averaged immediately for the sensor space analysis. For the source localization analysis, each session was taken through the source localization pipeline separately and then averaged in source space.

EEG sensor space analysis.

The average Fourier coefficients from each subject were averaged over a “sensor space region of interest” consisting of six electrodes over occipital cortex (70, 74, 75, 81, 82, 83). We used an additional second sensor space ROI, consisting of six electrodes over parietal cortex (53, 54, 61, 78, 79, 86) to generate control data. We isolated the response specific to symmetry by analyzing the odd and even harmonics of the spectrum separately (Norcia et al., 2002). Because we used a steady-state design in which P1 control images that had no rotation symmetry alternated with images from the four rotation symmetry groups, the rotation symmetry response will predominantly be found in the odd harmonics, whereas the even harmonics will mostly consist of responses arising from the image update, unrelated to symmetry. To demonstrate the specificity of the odd harmonics, we collected a separate dataset where we simply alternated the same P1 control images with an additional 10 P1 images (the same additional set of control images that was used in the fMRI control block). The expectation was that there will be no signal in the odd harmonics in this case, while the even harmonics will be largely unchanged (Norcia et al., 2002).

We computed the SNR for each harmonic based on the amplitude spectrum, the noise level for each harmonic, by taking the amplitude of one sideband on either side of the harmonic frequency over all four conditions, and averaging these eight values (two sidebands × four conditions) to get a single number representing the noise level across all conditions for that harmonic. We then computed SNR by dividing the amplitude at the harmonic frequency with the noise level. This was done individually for each participant. This procedure is analogous to the one used when computing SNR for the fMRI data, with the exception that we had to estimate the noise based on two, rather than four, side frequencies, to avoid including neighboring harmonics in the noise level.

Tissue segmentation procedure.

The FreeSurfer software package (http://surfer.nmr.mgh.harvard.edu) was used to extract both gray/white and gray/CSF boundaries, and generate cortical surface meshes. To avoid discontinuities in the cortically constrained inversion procedure arising from curvature differences between the gray/white and gray/CSF boundary, we generated a surface partway between these two boundaries that has gyri and sulci with approximately equal curvature. This “midgray” cortical surface consisted of a very dense triangular tessellation of several hundred thousand regularly spaced vertices. This tessellation was then downsampled to 20,484 vertices using the MNE software package (http://martinos.org/mne/stable/index.html). This number is low enough to compute the forward model on a standard workstation and yet accurately reflect the shape of cortical manifold (see e.g., Baillet et al., 2001). The final midgray surface was used to define the visual ROIs and the source space for the EEG current modeling. The FSL toolbox (http://www.fmrib.ox.ac.uk/fsl/) was used to segment, from the individual T1- and T2-weighted MRI scans, contiguous volume regions for the inner skull, outer skull, and scalp. These MRI volumes were then converted into inner skull, outer skull, and scalp surfaces (Smith, 2002; Smith et al., 2004) that define the boundaries between the brain/CSF and the skull, the skull and the scalp, and the scalp and the air.

EEG source-imaging analysis.

For each participant, the EEG source space was defined by the midgray surface (see the Tissue segmentation procedure section). The distance between the 20,484 vertices of this surface was on average 3.7 mm (SD, 1.5 mm; range, 0.1–11 mm). Current dipoles were placed at each of these vertices. Their orientations were constrained to be orthogonal to the cortical surface to diminish the number of parameters to be estimated in the inverse procedure (Hämäläinen et al., 1993). The source space, the 3D electrode locations, and the individually defined boundaries were then combined using the MNE software package to characterize the electric field propagation using a three-compartment boundary element method (Hämäläinen and Sarvas, 1989). The resulting forward model is linear and links the activity of the 20,484 cortical sources to the voltages recorded by our EEG electrodes.

Cortical current density estimates of the neural responses were obtained from an L2 minimum–norm inverse of the forward model as described in a recent study by Cottereau et al. (2012). We used the functionally defined visual ROIs to constrain these estimates by modifying the source–covariance matrix. The aim of this procedure was to decrease the tendency of the minimum–norm procedure to smooth activity across different functional ROIs. The following two modifications were applied: (1) we increased the variance allowed within the visual ROIs by a factor of two relative to other vertices; and (2) we enforced a local correlation constraint within each ROI using the first- and second-order neighborhoods on the cortical tessellation with a weighting function equal to 0.5 for the first order and 0.25 for the second order. This correlation constraint, therefore, respects both retinotopy and boundaries between visual areas, and permits a more precise dissociation of signals from different ROIs. This is not the case for other smoothing methods, such as LORETA, that apply the same smoothing rule across all of cortex (Pascual-Marqui et al., 1994). For each participant and condition, this inversion scheme was applied to the average Fourier coefficients (see the EEG acquisition and preprocessing section), which resulted in an estimation of these coefficients for each source on the cortical tessellation. The complex coefficients were then averaged across all the sources belonging to each functional ROI, for each participant. These within-participant ROI averages were then averaged across participants to generate cross-participant averages for each ROI. This procedure has been shown to reduce cross talk between different cortical ROIs, resulting in better area resolution in the group inversion than can be achieved in single participants (Cottereau et al., 2015).

Symmetry onset-timing estimation.

For the source localization analyses, our main interest was in estimating the response onset latency differences between different ROIs. We did this by fitting a piecewise linear function to the source-localized waveforms using a shape-language-modeling toolbox implemented in Matlab (http://www.mathworks.com/matlab-central/fileexchange/24443-slm-shape-language-modeling). The model was fit to the section of the waveform beginning at the onset of the second image (from one of the four wallpaper groups) and ending at 300 ms later, our rationale being that the response should begin within 300 ms of the image onset. The only constraint on the model was that it had to consist of three linear sections, separated by two breakpoints that were determined by the fitting procedure. The first breakpoint provided a reasonable estimate of the latency of the negative deflection characteristic of the symmetry response (see Fig. 9).

We used a jackknife procedure to estimate the error of the latency differences between pairs of ROIs by computing average waveforms across all participants except one, and then fitting the piecewise linear model to the resulting waveform. We then subtracted the resulting response latency in the first ROI from the latency in the second. We repeated this procedure for all participants, and estimated the SE of the differences, as described by Miller et al. (1998, their Eq. 2).

Computational modeling.

Wallpaper patterns can be construed as being textures or, alternatively, as having elements of form that relate to object properties. There has been considerable interest in developing computational models of early visual cortical mechanisms or machine vision algorithms that could support either texture or object processing. When interpreting our findings, it is important to assess the extent to which differences in the neural responses elicited by the four wallpaper groups can be explained by properties of the stimuli that are unrelated to rotation symmetry, per se, but could be relatable to other types of regularity found in textures. To address this question, we measured the predicted strength of neural responses to the four wallpaper groups according to (1) their Weibull statistics, (2) their discriminability by a robust texture classifier, (3) their pattern phase congruency across spatial scale, and the (4) the responses of the model of second-order contrast sensitivity. Each of these approaches has been shown to be sensitive to the statistics of natural images and can provide quantitative metrics of texture similarity. Our computational approach was analogous to the one taken in the EEG and fMRI experiments, in that we assessed the predicted response difference between each group of exemplars and its corresponding P1 control. All four models were run on the same set of images that was used for the EEG and fMRI experiments, with 20 exemplars per wallpaper group, and the predicted responses differences were then averaged across exemplars within a group. The first three approaches were applied to cropped versions of the stimuli to reduce the number of artifacts caused by the aperture, while the contrast model was applied to the whole image, including the aperture, and the resulting output map was then cropped.

It has been suggested that parameters of the Weibull distribution provide useful summary statistics for natural images (Geusebroek and Smeulders, 2002). The Weibull distribution was initially used to describe the range of particle sizes resulting from roller mills (Rosin and Rammler, 1933), and it has intuitive appeal as an index of “fragmentation” in an image (Geusebroek and Smeulders, 2002; Groen et al., 2013). The Weibull distribution is parameterized by a “scale” parameter a and a “shape” variable b. The scale parameter varies with the contrast energy of the image, and the shape parameter varies with the spatial coherence of the image (Groen et al., 2013). Several recent event-related potential (ERP) reports have investigated the extent to which Weibull statistics could explain the variance in the neural responses (Scholte et al., 2009; Groen et al., 2012a,b; Groen et al., 2013). The Weibull statistics capture a substantial fraction of the ERP variance in response to natural images that vary widely in their structure. We computed the Weibull statistics (a and b) for each image, I(x,y), as follows: First, the image was convolved with first-order Gaussian derivative filters for some Gaussian G(x, y; σ). We then computed the gradient magnitude over the image as follows: Then, a 256-bin histogram was calculated over ∇¹I, which was modeled using a Weibull distribution, as follows: using the wblfit command in MATLAB. To maximize the fit between the Weibull distribution and our data, a and b were calculated for a range of spatial scales (2 ≤ σ ≤ 6), and we used the scale that gave the highest (logistic regression) accuracy when classifying the exemplars from their matched P1 controls. Finally, we used β_a and β_b to scale a and b and construct a distance metric between pairs of images I₁ and I₂, as follows: where a_i, b_i are the Weibull statistics from image I_i. This distance metric was generated for each exemplar and its matched P1 control.

Our second approach was to ask whether features identified as informative by the computer vision and texture classification literature can predict the neural responses to our symmetry patterns. A recent study (Dong et al., 2014) evaluated 51 computational texture features in terms of how well they could rank patterns in terms of their visual similarity and found that classification based on local image-patch features (Varma and Zisserman, 2003) offered the best performance. The local image-patch approach is also attractive because it shares several similarities with the sparse coding model of Olshausen and Field (1997). In particular, both models are based on extracting n × n pixel local neighborhoods (“textons”) from a set of images and then using an unsupervised clustering algorithm to learn a set of textons that can be used to represent the images. The main difference between approaches is the algorithm used to learn the texton dictionary, although both methods learn dictionaries of a similar size. Varma and Zisserman (2003) used dictionaries consisting of 1000–3000 exemplar textons, while Olshausen (2013) has recently experimented with 10-fold overcomplete sparse coding dictionaries consisting of 2048 basis functions. Critically, as noted above, features from the local patch model were among the best predictors of human perceptual similarity ratings for a large texture database (Clarke et al., 2011a). Our implementation involves downsampling the image by a factor of four, and then using k-means to learn a dictionary of (k = 100) n × n local neighborhoods (n = 3). Once a texton dictionary is created, an image, I_i, can be represented by a histogram, h_i, over the textons, and two textures are compared with one another using a simplified form of the Bhattacharyya distance, as follows: Using this approach, we calculated the distance between each exemplar and its corresponding P1 control.

Natural images contain numerous high-order correlations, among which correlations in the phase spectrum across scale are particularly prominent. Several models of human feature perception have been built to extract features in the frequency domain that are due to the presence of various types of edges that can create phase correlations. Inherent in these models is the extraction of correlations in spatial phase across scale (Morrone and Burr, 1988; Owens, 1994; Kovesi, 2000). Perceptual appearance is strongly controlled by cross-scale phase coherence (Morrone and Burr, 1988), and sensitivity to certain changes in the appearance of a texture can be accurately modeled (R² > 0.95) using the variance of the phase congruency map (Emrith et al., 2010). We tested whether patterns of phase coherency across scale could predict the differences in neural responses between the wallpaper groups, by computing a summary difference on phase coherency in terms of the log of the variance of the phase congruency map between each exemplar image and its matched P1 control.

The second-order contrast (SOC) model is a cascaded, feedforward model of BOLD responses to contrast and texture, in which the image is first processed through a bank of contrast-normalized, localized, V1-like filters, the output of which is then reprocessed by a second stage that, among other things, measures the contrast variability in the output of the first stage (Kay et al., 2013). The SOC model has been shown to be effective at predicting differences in texture responses between V1 and V2/3 (Kay et al., 2013), so it is a plausible candidate model for predicting differences in the response to the different wallpaper groups. Furthermore, the model is attractive as a test case because second-order contrast is higher in natural images with intact phase spectra (Kay et al., 2013), and because the model has been fit to data for V1, V2, V3, and V4, successfully capturing increased sensitivity to the structure of natural scenes in the higher-order areas. We applied the SOC model separately to each exemplar image, using model parameters for V1, V2, V3, and V4 provided in the original publication (Kay et al., 2013). We then summed the model output across each image, and computed the modeled difference between each exemplar and its matched P1 control.