diff --git a/perlstim/_noise.h b/perlstim/_noise.h deleted file mode 100644 index 2a97d65..0000000 --- a/perlstim/_noise.h +++ /dev/null @@ -1,107 +0,0 @@ -// Copyright (c) 2008, Casey Duncan (casey dot duncan at gmail dot com) -// see LICENSE.txt for details - -#include - -#if defined(_MSC_VER) -#define inline _inline -#define M_1_PI 0.31830988618379067154 -#endif - -const float GRAD3[][3] = { - {1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0}, - {1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1}, - {0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}, - {1,0,-1},{-1,0,-1},{0,-1,1},{0,1,1}}; - -const float GRAD4[][4] = { - {0,1,1,1}, {0,1,1,-1}, {0,1,-1,1}, {0,1,-1,-1}, - {0,-1,1,1}, {0,-1,1,-1}, {0,-1,-1,1}, {0,-1,-1,-1}, - {1,0,1,1}, {1,0,1,-1}, {1,0,-1,1}, {1,0,-1,-1}, - {-1,0,1,1}, {-1,0,1,-1}, {-1,0,-1,1}, {-1,0,-1,-1}, - {1,1,0,1}, {1,1,0,-1}, {1,-1,0,1}, {1,-1,0,-1}, - {-1,1,0,1}, {-1,1,0,-1}, {-1,-1,0,1}, {-1,-1,0,-1}, - {1,1,1,0}, {1,1,-1,0}, {1,-1,1,0}, {1,-1,-1,0}, - {-1,1,1,0}, {-1,1,-1,0}, {-1,-1,1,0}, {-1,-1,-1,0}}; - -// At the possible cost of unaligned access, we use char instead of -// int here to try to ensure that this table fits in L1 cache -const unsigned char PERM[] = { - 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, - 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, - 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, - 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, - 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, - 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, - 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, - 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, - 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, - 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, - 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, - 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, - 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, - 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, - 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, - 141, 128, 195, 78, 66, 215, 61, 156, 180, 151, 160, 137, 91, 90, 15, 131, - 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, - 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, - 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, - 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, - 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, - 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, - 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, - 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, - 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, - 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, - 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, - 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, - 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, - 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, - 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, - 180}; - -const unsigned char SIMPLEX[][4] = { - {0,1,2,3},{0,1,3,2},{0,0,0,0},{0,2,3,1},{0,0,0,0},{0,0,0,0},{0,0,0,0}, - {1,2,3,0},{0,2,1,3},{0,0,0,0},{0,3,1,2},{0,3,2,1},{0,0,0,0},{0,0,0,0}, - {0,0,0,0},{1,3,2,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}, - {0,0,0,0},{0,0,0,0},{0,0,0,0},{1,2,0,3},{0,0,0,0},{1,3,0,2},{0,0,0,0}, - {0,0,0,0},{0,0,0,0},{2,3,0,1},{2,3,1,0},{1,0,2,3},{1,0,3,2},{0,0,0,0}, - {0,0,0,0},{0,0,0,0},{2,0,3,1},{0,0,0,0},{2,1,3,0},{0,0,0,0},{0,0,0,0}, - {0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0},{2,0,1,3}, - {0,0,0,0},{0,0,0,0},{0,0,0,0},{3,0,1,2},{3,0,2,1},{0,0,0,0},{3,1,2,0}, - {2,1,0,3},{0,0,0,0},{0,0,0,0},{0,0,0,0},{3,1,0,2},{0,0,0,0},{3,2,0,1}, - {3,2,1,0}}; - -#define fastfloor(n) (int)(n) - (((n) < 0.0f) & ((n) != (int)(n))) - -// Fast sine/cosine functions from -// http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/page__st__80 -// Note the input to these functions is not radians -// instead x = [0, 2] for r = [0, 2*PI] - -static inline float fast_sin(float x) -{ - // Convert the input value to a range of -1 to 1 - // x = x * (1.0f / PI); - - // Wrap around - volatile float z = (x + 25165824.0f); - x = x - (z - 25165824.0f); - - #if LOW_SINE_PRECISION - return 4.0f * (x - x * fabsf(x)); - #else - { - float y = x - x * fabsf(x); - const float Q = 3.1f; - const float P = 3.6f; - return y * (Q + P * fabsf(y)); - } - #endif -} - -static inline float fast_cos(float x) -{ - return fast_sin(x + 0.5f); -} - diff --git a/setup.py b/setup.py index d4dad6d..3bee3db 100644 --- a/setup.py +++ b/setup.py @@ -8,27 +8,28 @@ if sys.platform != 'win32': compile_args = ['-funroll-loops'] -else: - # XXX insert win32 flag to unroll loops here - compile_args = [] + +with open("README.md", "r") as f: + long_desc = f.read() setup( - name='perlstim', + name='zebranoise', version='0.1', - description='Perlin noise stimuli for Python', - long_description='Based on the "noise" package by Casey Duncan', + description='Zebra noise stimuli for Python', + long_description = long_desc, + long_description_content_type='text/markdown', author='Max Shinn', author_email='m.shinn@ucl.ac.uk', - url='https://github.com/mwshinn/perlin_stimuli', + url='https://github.com/mwshinn/zebra_noise', classifiers = [ 'Development Status :: 4 - Beta', 'Programming Language :: C', 'Programming Language :: Python :: 3', ], - - packages=['perlstim'], + install_requires = ['numpy', 'imageio', 'imageio-ffmpeg'], + packages=['zebranoise'], ext_modules=[ - Extension('perlstim._perlin', ['perlstim/_perlin.c'], + Extension('zebranoise._perlin', ['zebranoise/_perlin.c'], extra_compile_args=compile_args, include_dirs=[np.get_include()], ) diff --git a/zebranoise/__init__.py b/zebranoise/__init__.py new file mode 100644 index 0000000..72161b7 --- /dev/null +++ b/zebranoise/__init__.py @@ -0,0 +1,4 @@ +from .perlin_stimulus import PerlinStimulus +from .util import generate_frames + +Perl = PerlinStimulus # Backward compatibility diff --git a/perlstim/_perlin.c b/zebranoise/_perlin.c similarity index 68% rename from perlstim/_perlin.c rename to zebranoise/_perlin.c index 80baac4..0b47cf2 100644 --- a/perlstim/_perlin.c +++ b/zebranoise/_perlin.c @@ -5,7 +5,6 @@ #include "Python.h" #include #include -#include "_noise.h" #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include @@ -15,6 +14,49 @@ #define lerp(t, a, b) ((a) + (t) * ((b) - (a))) +const float GRAD3[][3] = { + {1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0}, + {1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1}, + {0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}, + {1,0,-1},{-1,0,-1},{0,-1,1},{0,1,1}}; + +// At the possible cost of unaligned access, we use char instead of +// int here to try to ensure that this table fits in L1 cache +const unsigned char PERM[] = { + 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, + 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, + 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, + 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, + 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, + 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, + 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, + 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, + 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, + 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, + 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, + 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, + 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, + 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, + 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, + 141, 128, 195, 78, 66, 215, 61, 156, 180, 151, 160, 137, 91, 90, 15, 131, + 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, + 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, + 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, + 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, + 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, + 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, + 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, + 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, + 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, + 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, + 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, + 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, + 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, + 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, + 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, + 180}; + + static inline float grad3(const int hash, const float x, const float y, const float z) { diff --git a/perlstim/__init__.py b/zebranoise/perlin_stimulus.py similarity index 61% rename from perlstim/__init__.py rename to zebranoise/perlin_stimulus.py index 319b3ec..50b2259 100644 --- a/perlstim/__init__.py +++ b/zebranoise/perlin_stimulus.py @@ -1,22 +1,30 @@ -from subprocess import call -import imageio import hashlib import os -from . import _perlin -import numpy as np import warnings -import scipy.ndimage -from pathlib import Path import tempfile +from subprocess import call +from pathlib import Path +import numpy as np +import imageio +from imageio_ffmpeg import get_ffmpeg_exe -class Perl: +from . import _perlin +from .util import generate_frames, filter_frames, filter_frames_index_function, XYSCALEBASE + + +class PerlinStimulus: """Class to generate a Perlin noise stimulus. + This is the "heavy duty" way of generating stimuli, mostly for experimenting + with new Perlin-based stimuli. If you would just like to generate zebra + noise, use the "zebra_noise" function instead, which is faster and more + memory efficient. + To generate noise, instantiate this method. This will automatically generate Perlin noise which you can save, filter, etc. + """ - XYSCALEBASE = 100 - def __init__(self, xsize, ysize, tdur, levels=10, xyscale=.5, tscale=1, fps=30, xscale=1.0, yscale=1.0, seed=0, demean="both", cachedir="perlcache", delay_batch=False): + def __init__(self, xsize, ysize, tdur, levels=10, xyscale=.2, tscale=50, fps=30, xscale=1.0, yscale=1.0, seed=0, demean="both", cachedir="perlcache", delay_batch=False): """Initialise Perlin noise stimulus. Parameters @@ -54,8 +62,8 @@ def __init__(self, xsize, ysize, tdur, levels=10, xyscale=.5, tscale=1, fps=30, """ assert demean in ["both", "time", "space", "none"] tsize = int(tdur*fps) - self.ratio = xsize/ysize*self.XYSCALEBASE - assert self.ratio == int(self.ratio), "Ratio between x and y times 100 must be an integer" + self.ratio = xsize/ysize*XYSCALEBASE + assert self.ratio == int(self.ratio), f"Ratio between x and y times {XYSCALEBASE} must be an integer" self.ratio = int(self.ratio) textra = (tscale - (tsize % tscale)) % tscale if textra > 0: @@ -118,104 +126,6 @@ def discretize(cls, im): im *= 255 ret = im.astype(np.uint8) return ret - @classmethod - def filter(cls, im, filt, *args): - """Apply a filter to a batch - - Parameters - ---------- - im : 3D float ndarray, values ∈ [0,1] - Noise movie - filt : str - The name of the filter - *args : tuple - Extra arguments are passed to the filter - - Returns - ------- - im : 3D float ndarray, values ∈ [0,1] - Filtered noise movie - """ - if filt == "threshold": - return (im>args[0]).astype(np.float32) - if filt == "softthresh": - return 1/(1+np.exp(-args[0]*(im-.5))) - if filt == "comb": - return (im//args[0] % 2 == 1).astype(np.float32) - if filt == "invert": - return 1-im - if filt == "reverse": - return im # We need to use filter_index_function for this - if filt == "blur": - return np.asarray([scipy.ndimage.filters.gaussian_filter(im.astype(np.float32)[:,:,i], args[0], mode='wrap') for i in range(0, im.shape[2])]).transpose([1,2,0]) - if filt == "wood": - return (im % args[0]) / args[0] - if filt == "center": - return 1-(np.abs(im-.5)*2) - if filt == "photodiode": - im = im.copy() - s = args[0] - im[:s,-s:,::2] = 0 - im[:s,-s:,1::2] = 1 - return im - if filt == "photodiode_anywhere": - im = im.copy() - x = args[0] - y = args[1] - s = args[2] - im[y:(y+s),x:(x+s),::2] = 0 - im[y:(y+s),x:(x+s),1::2] = 1 - return im - if filt == "photodiode_b2": - im = im.copy() - s = 125 - im[:s,-s:,::2] = 0 - im[:s,-s:,1::2] = 1 - return im - if filt == "photodiode_fusi": - im = im.copy() - s = 75 - im[:s,-s:,::2] = 0 - im[:s,-s:,1::2] = 1 - return im - if filt == "photodiode_bscope": - im = im.copy() - s = 100 - im[-s:,:s,::2] = 0 - im[-s:,:s,1::2] = 1 - return im - if callable(filt): - return filt(im) - raise ValueError("Invalid filter specified") - def filter_index_function(self, filters): - """Reordering frames in the video based on the filter. - - - Parameters - ---------- - filters : list of strings or tuples - the list of filters passed to save_video - - Returns - ------- - function mapping int -> int - Reindexing function - - Notes - ----- - Some filters may need to operate on the global video instead of in - batches. However, for large videos, batches are necessary due to - limited amounts of RAM. Thus, this function should return another - function which takes an index as input and outputs a new index, - remapping the initial noise frame to the output video frame. This - was primarily designed to support reversing the video, but it might be - useful for other things too. - - """ - if "reverse" in filters: - return lambda x : self.nframes - x - 1 - return lambda x : x - def generate_frame(self, t=0, filters=[]): """Generate and return a single image of noise. @@ -223,7 +133,7 @@ def generate_frame(self, t=0, filters=[]): Parameters ---------- t : float - the timepoint at which to generate the image + the frame at which to generate the image filters : list of str and/or (str, ...) tuples A list of filters to apply to the stimulus. If a filter requires parameters, pass a tuple, where the first element is the name of @@ -234,17 +144,9 @@ def generate_frame(self, t=0, filters=[]): 2-dimensional ndarray An image of the noise """ - tunits = int(self.size[2]/self.tscale) - arr = _perlin.make_perlin(np.arange(0, self.size[0], dtype="float32")/self.size[1]/self.xscale, # Yes, divide by y size - np.arange(0, self.size[1], dtype="float32")/self.size[1]/self.yscale, - np.asarray([t], dtype="float32"), - octaves=self.levels, - persistence=self.xyscale, - repeatx=self.ratio, - repeaty=self.XYSCALEBASE, - repeatz=tunits, - base=self.seed) - arr = arr.swapaxes(0,1) + arr = generate_frames(self.size[0], self.size[1], self.size[2], timepoints=[t], levels=self.levels, + xyscale=self.xyscale, tscale=self.tscale, xscale=self.xscale, + yscale=self.yscale, seed=self.seed) if self.demean in ["both", "time"]: arr -= np.mean(arr, axis=(0,1), keepdims=True) for f in filters: @@ -254,7 +156,7 @@ def generate_frame(self, t=0, filters=[]): else: n = f[0] args = f[1:] - arr = self.filter(arr, n, *args) + arr = filter_frames(arr, n, *args) return arr.squeeze() def generate_batch(self): @@ -274,10 +176,6 @@ def generate_batch(self): self.max_ = stats['max_'] self.nframes = stats['nframes'] return - # Use the temporal scale and number of timepoints to compute how many - # units to make the stimulus across the temporal dimension - tunits = int(self.size[2]/self.tscale) - ts_all = np.arange(0, self.size[2], dtype="float32")/self.tscale # Generate the stimuli and save the means and mins/maxes. If we are # demeaning across space, we won't end up using the mins or maxes saved # here. @@ -286,19 +184,12 @@ def generate_batch(self): weights = [] mins = [] maxes = [] - for k in range(0, int(np.ceil(len(ts_all)/self.batch_size))): - ts = ts_all[(k*self.batch_size):((k+1)*self.batch_size)] - print("batch", k, len(ts_all), self.batch_size, len(ts)) - arr = _perlin.make_perlin(np.arange(0, self.size[0], dtype="float32")/self.size[1]/self.xscale, # Yes, divide by y size - np.arange(0, self.size[1], dtype="float32")/self.size[1]/self.yscale, - ts, - octaves=self.levels, - persistence=self.xyscale, - repeatx=self.ratio, - repeaty=self.XYSCALEBASE, - repeatz=tunits, - base=self.seed) - arr = arr.swapaxes(0,1) + for k in range(0, int(np.ceil(self.size[2]/self.batch_size))): + print(f"Generating batch {k}") + arr = generate_frames(self.size[0], self.size[1], tsize=self.size[2], + timepoints=list(range(k*self.batch_size,min(self.size[2],(k+1)*self.batch_size))), + levels=self.levels, xyscale=self.xyscale, tscale=self.tscale, xscale=self.xscale, + yscale=self.yscale, seed=self.seed) if self.demean in ["both", "time"]: arr -= np.mean(arr, axis=(0,1), keepdims=True) mins.append(np.min(arr)) @@ -340,9 +231,10 @@ def save_grey_pad(self, fn, dur, bitrate=20): os.link(self.tmpdir.joinpath("_grey.tif"), self.tmpdir.joinpath(f"_frame{i:05}.tif")) if fn[-4:] != ".mp4": fn += ".mp4" - call(["ffmpeg", "-r", str(self.fps), "-i", str(self.tmpdir.joinpath("_frame%5d.tif")), "-c:v", "mpeg2video", "-an", "-b:v", f"{bitrate}M", fn]) - call([f"rm "+str(self.tmpdir.joinpath("_frame*.tif"))], shell=True) - call([f"rm "+str(self.tmpdir.joinpath("_grey.tif"))], shell=True) + call([get_ffmpeg_exe(), "-r", str(self.fps), "-i", str(self.tmpdir.joinpath("_frame%5d.tif")), "-c:v", "mpeg2video", "-an", "-b:v", f"{bitrate}M", fn]) + for p in self.tmpdir.glob("_frame*.tif"): + p.unlink() + self.tmpdir.joinpath("_grey.tif").unlink() def save_video(self, fn, loop=1, filters=[], bitrate=20): """Save the filtered stimulus. @@ -365,7 +257,7 @@ def save_video(self, fn, loop=1, filters=[], bitrate=20): raise IOError("Output video file already exists!") i = 0 k = 0 - filtind = self.filter_index_function(filters) + filtind = filter_frames_index_function(filters, nframes=self.nframes) while os.path.isfile(self.cache_filename(k)): data = np.load(self.cache_filename(k)).astype("float32") # Renormalise using precomputed mins/maxes @@ -379,7 +271,7 @@ def save_video(self, fn, loop=1, filters=[], bitrate=20): else: n = f[0] args = f[1:] - data = self.filter(data, n, *args) + data = filter_frames(data, n, *args) data = self.discretize(data) assert data.dtype == 'uint8' for j in range(0, data.shape[2]): @@ -393,5 +285,6 @@ def save_video(self, fn, loop=1, filters=[], bitrate=20): os.link(self.tmpdir.joinpath(f"_frame{i:05}.tif"), self.tmpdir.joinpath(f"_frame{(i+j*n_frames):05}.tif")) if fn[-4:] != ".mp4": fn += ".mp4" - call(["ffmpeg", "-r", str(self.fps), "-i", str(self.tmpdir.joinpath("_frame%5d.tif")), "-c:v", "mpeg2video", "-an", "-b:v", f"{bitrate}M", fn]) - call([f"rm "+str(self.tmpdir.joinpath("_frame*.tif"))], shell=True) + call([get_ffmpeg_exe(), "-r", str(self.fps), "-i", str(self.tmpdir.joinpath("_frame%5d.tif")), "-c:v", "mpeg2video", "-an", "-b:v", f"{bitrate}M", fn]) + for p in self.tmpdir.glob("_frame*.tif"): + p.unlink() diff --git a/zebranoise/util.py b/zebranoise/util.py new file mode 100644 index 0000000..fac770d --- /dev/null +++ b/zebranoise/util.py @@ -0,0 +1,123 @@ +import numpy as np + +from . import _perlin + +XYSCALEBASE = 100 + +def filter_frames(im, filt, *args): + """Apply a filter/transformation to an image batch + + Parameters + ---------- + im : 3D float ndarray, values ∈ [0,1] + Frames to filter + filt : str + The name of the filter + *args : tuple + Extra arguments are passed to the filter + + Returns + ------- + im : 3D float ndarray, values ∈ [0,1] + Filtered noise movie + """ + if filt == "threshold": + return (im>args[0]).astype(np.float32) + if filt == "softthresh": + return 1/(1+np.exp(-args[0]*(im-.5))) + if filt == "comb": + return (im//args[0] % 2 == 1).astype(np.float32) + if filt == "invert": + return 1-im + if filt == "reverse": + return im # We need to use filter_index_function for this + if filt == "blur": + return np.asarray([scipy.ndimage.filters.gaussian_filter(im.astype(np.float32)[:,:,i], args[0], mode='wrap') for i in range(0, im.shape[2])]).transpose([1,2,0]) + if filt == "wood": + return (im % args[0]) / args[0] + if filt == "center": + return 1-(np.abs(im-.5)*2) + if filt == "photodiode": + im = im.copy() + s = args[0] + im[:s,-s:,::2] = 0 + im[:s,-s:,1::2] = 1 + return im + if filt == "photodiode_anywhere": + im = im.copy() + x = args[0] + y = args[1] + s = args[2] + im[y:(y+s),x:(x+s),::2] = 0 + im[y:(y+s),x:(x+s),1::2] = 1 + return im + if filt == "photodiode_b2": + im = im.copy() + s = 125 + im[:s,-s:,::2] = 0 + im[:s,-s:,1::2] = 1 + return im + if filt == "photodiode_fusi": + im = im.copy() + s = 75 + im[:s,-s:,::2] = 0 + im[:s,-s:,1::2] = 1 + return im + if filt == "photodiode_bscope": + im = im.copy() + s = 100 + im[-s:,:s,::2] = 0 + im[-s:,:s,1::2] = 1 + return im + if callable(filt): + return filt(im) + raise ValueError("Invalid filter specified") + +def filter_frames_index_function(filters, nframes): + """Reordering frames in the video based on the filter. + + + Parameters + ---------- + filters : list of strings or tuples + the list of filters passed to save_video + + Returns + ------- + function mapping int -> int + Reindexing function + + Notes + ----- + Some filters may need to operate on the global video instead of in + batches. However, for large videos, batches are necessary due to + limited amounts of RAM. Thus, this function should return another + function which takes an index as input and outputs a new index, + remapping the initial noise frame to the output video frame. This + was primarily designed to support reversing the video, but it might be + useful for other things too. + + """ + if "reverse" in filters: + return lambda x : nframes - x - 1 + return lambda x : x + +def generate_frames(xsize, ysize, tsize, timepoints, levels=10, xyscale=.5, tscale=1, xscale=1.0, yscale=1.0, seed=0): + """Preprocess arguments before passing to the C implementation of Perlin noise. + """ + # Use the temporal scale and number of timepoints to compute how many + # units to make the stimulus across the temporal dimension + tunits = int(tsize/tscale) + ts_all = np.arange(0, tsize, dtype="float32")/tscale + ratio = int(xsize/ysize*XYSCALEBASE) + arr = _perlin.make_perlin(np.arange(0, xsize, dtype="float32")/ysize/xscale, # Yes, divide by y size + np.arange(0, ysize, dtype="float32")/ysize/yscale, + ts_all[timepoints], + octaves=levels, + persistence=xyscale, + repeatx=ratio, + repeaty=XYSCALEBASE, + repeatz=tunits, + base=seed) + arr = arr.swapaxes(0,1) + return arr