# A domain adaptive deep learning solution for scanpath prediction of
paintings

Mohamed Amine Kerkouri*<EMAIL_ADDRESS>Université
d’Orléans12 Rue de BloisOrléansFrance45100 ,  Marouane Tliba
<EMAIL_ADDRESS>Université d’Orléans12 Rue de
BloisOrléansFrance45100 ,  Aladine Chetouani  Université d’Orléans12 Rue de
BloisOrléansFrance45100 and  Alessandro Bruno  Bournemouth UniversityTalbot
Campus, Fern BarrowPooleUnited KingdomBH125BB

###### Abstract.

Cultural heritage understanding and preservation is an important issue for
society as it represents a fundamental aspect of its identity. Paintings
represent a significant part of cultural heritage, and are the subject of
study continuously. However, the way viewers perceive paintings is strictly
related to the so-called HVS (Human Vision System) behaviour. This paper
focuses on the eye-movement analysis of viewers during the visual experience
of a certain number of paintings. In further details, we introduce a new
approach to predicting human visual attention, which impacts several cognitive
functions for humans, including the fundamental understanding of a scene, and
then extend it to painting images. The proposed new architecture ingests
images and returns scanpaths, a sequence of points featuring a high likelihood
of catching viewers’ attention. We use an FCNN (Fully Convolutional Neural
Network), in which we exploit a differentiable channel-wise selection and
Soft-Argmax modules. We also incorporate learnable Gaussian distributions onto
the network bottleneck to simulate visual attention process bias in natural
scene images. Furthermore, to reduce the effect of shifts between different
domains (i.e. natural images, painting), we urge the model to learn
unsupervised general features from other domains using a gradient reversal
classifier. The results obtained by our model outperform existing state-of-
the-art ones in terms of accuracy and efficiency.

Scanpath Prediction, Unsupervised Domain Adaptation, Paintings.

* Equal contribution; the order of first authors was randomly selected

##  1\. Introduction

Cultural heritage encompasses the habits, traditions, artifacts, and artistic
expressions produced by a community and carried down through generations
(Hartmann, 2020). Thus representing an integral part of their communal
identity which should be preserved.

Paintings represent an important part of the artistic heritage, it is a
picture created by putting paint on a surface, or the activity or skill of
creating pictures by using
paint111https://dictionary.cambridge.org/dictionary/english/painting. They
depict a pictorial recounting of events, ideas, and concepts of the
environment, era, and tendencies of the artist. Thus, the preservation and
understanding of such heritage are important as they are an integral part of
human history.

Understanding this type of media is directly linked to the visual perception
mechanisms of the human visual system (HVS) and as human eyes have a
resolution of more than 500 MP and receive about 10 Gigabits of information
per second (Raichle, 2010).

Due to the enormous quantity of data to be processed, HVS needs to filter the
most relevant areas in a visual scene while discarding the least meaningful
ones from a perceptual perspective.

In greater detail, eye movements allow filtering out non-relevant information
by attending to a region of the visual field rather than others. That
mechanism unburdens the cognitive load and increases the efficiency of the
visual perception system. Findings revealed that low-level visual features
(color, edges, texture, intensity, contrast, etc) catch eye movements in the
early stages of image observation. In contrast, semantic visual features
(faces, text, humans, objects, etc.) play a critical role in a later stage of
the visual experience (Treisman and Gelade, 1980).

The aforementioned mechanism can be categorized into ’covert’ and ’overt’
attention. The former focuses on optimizing the processing activity in the
brain cortex. The latter has a strict relation with eye movements that bring
regions of interest of the visual field onto the fovea. Eye-trackers capture
eye movements and allow drawing the gaze trajectory of viewers while exploring
a specific visual stimulus. The scanpath consists of fixation points and
saccades. Saccades are fast eye movements with no visual perception between
fixation points (Goldberg and Helfman, 2010).

Fixations of multiple observers for a specific stimulus (an image or a video
sequence) are generally represented with the so-called fixation point maps.
These maps can be smoothed using a Gaussian kernel representing $1^{\circ}$ of
visual angle to give a spatial distribution of eye movements. Furthermore,
they provide a probability heat map called a saliency map where each pixel has
a normalized probability value of grabbing the viewer’s attention.

The main goal of saliency models is to generate saliency maps as close as
possible to fixation point maps (Tliba et al., 2022a), (Tliba et al., 2022b).

Attention can also be classified using the ”neural processing path” it uses
for the mechanism. ”Bottom-up Attention” is mainly directed by the external
stimuli features, like the one aforementioned previously. This means that the
intrinsic properties and features of the images reign over the mechanism, and
this correlates with the theory that states that attention is pre-processing
operation. This type is unintentional, fast and instinctual.

Meanwhile, the other type is called ”Top-down Attention”. This type is
influenced by other highly cognitive processes originating from the prefrontal
cortex and the long-term memory, it relates to the theory stating that
attention is a post-processing operation. This mechanism is intentional, slow,
and premeditated.

The output of saliency models is used in several topics such as image quality
assessment (Chetouani, 2018b, a; Abouelaziz et al., 2019, 2020; Chetouani and
Li, 2020), image and video compression (Patel et al., 2021), image captioning
and description (Cornia et al., 2017), image search and retrieval (Wang et
al., 2020), image enhancement for people with CVD (Colour Vision Deficiency)
(Bruno et al., 2019), saliency led painting restoration (Yue and Jun, 2021)
and so forth (et al, 2019).

Scientific research in saliency prediction was pioneered by the seminal work
of (Koch and Ullman, 1987), implemented by (Itti and Koch, 2001) using a
multi-scale model reliant on the extraction of low-level features (i.e. color,
intensity, and orientation).

Alongside visual saliency, scanpath prediction has gained the interest of
researchers lately. The winner-take-all (WTA) module presented by (Itti and
Koch, 2001) was the first to predict scanpaths. They extracted the points with
the highest saliency value as fixation points and inhibited the region
surrounding that point. In (Le Meur and Liu, 2015), the authors proposed a
stochastic approach to generate scanpaths. They predicted saliency maps and
modeled the probabilities of several biases (i.e. saccade amplitudes and
saccade orientations).

The authors of (Zanca et al., 2019) analogized the saliency map with a 2D
gravitational field affecting the trajectory of a mass representing the gaze.
The authors in (Assens et al., 2018) developed a Long-Short-Term-Memory (LSTM)
network. They combined it with a conditional GAN (Generative Adversarial
Network) architecture to predict the scanpath of a given visual stimulus. A
framework employing foveal saliency maps, temporal duration, and a model of
the Inhibition of Return (IoR) to predict scanpath was proposed in (Bao and
Chen, 2020).

The authors of (Kerkouri et al., 2021) introduced an end-to-end model which
simultaneously predicts the scanpath and the saliency map of an image. Then,
in (Kerkouri and Chetouani, 2021) they used a fully convolutional architecture
of known computer vision models in a regressive manner to predict scanpaths.

Figure 1. General architecture of the proposed method.

In this paper, we propose a novel efficient method for vector to sequence
modeling (Image to scanpath) using fully Convolutional Neural Networks.
Because images originating from different domains ( i.e. natural scenes, art,
synthetic, … etc) exhibit different biases related to the domains themselves
and at the same time visual perception tasks also have their own subjective
biases. It becomes hard for models to model relevant features specific to the
images while ignoring these types of noisy distributions. To that end, We
investigate two sorts of domain adaptation techniques one related to the
downstream task priors biases formulation and the other focusing on feature
extraction generalization also known as domain representation alignment, over
different data distribution domains and employ this technique to improve the
prediction on paintings dataset.

The main contributions of this work are summarized as it follows:

  * •

Propose a novel and efficient deep model for scanpath prediction.

  * •

Investigate domain adaptation on scanpath prediction by learning human visual
attention bias and use unsupervised adversarial domain adaption to generalize
our model prediction on other image domains ( paintings, etc.).

  * •

Conduct an extensive qualitative and quantitative experimental campaign to the
proposed method’s effectiveness.

The rest of the paper is organized as follows: in Section 2, we describe the
components of the proposed model architecture as well as the training protocol
and the generalization process using the unsupervised domain adaptation. In
section 3, the results and the evaluation of the model are detailed. Section 4
ends the paper with conclusions.

##  2\. Proposed Method

This work aims to introduce a new framework for scanpath prediction and
demonstrate the effectiveness of unsupervised domain adaptation on the task.
The general architecture of the proposed model is graphically depicted in
Fig.1. We first use MobileNet as the backbone for feature extraction.
Secondly, learnable domain prior bias maps are concatenated to the output of
the feature extractor and then passed onto a fully convolutional decoder. The
latter is composed of several modules that aim to extract specific
characteristics relative to the final prediction. Global pooling features are
extracted from the maps resulting from the merging module and then passed to a
fully connected network in order to predict a probability vector. The latter
is then discretized into a binary mask related to the relevant predicted
fixation points of the scanpath. The probability vector is used as a channel-
wise weighting vector to accentuate the importance of specific channels
compared to others. Finally, the resulting weighted feature maps are passed
through a Soft-ArgMax function, where each fixation point coordinate is
predicted from each feature map.222code will soon be availble at :
https://github.com/Submit-code/ScanpathDomain

###  2.1. Feature extractor

Pixel-wise computations can be labor-intensive on dimensional complexity and
spatial correlation. CNNs are used to model their features before extracting
them. Different architectures were proposed and evolved for this end (Mahmoudi
et al., 2020a, b). Here, we use MobileNet (Sandler et al., 2018) to extract
features from the input images. It is a lightweight CNN built for use in
mobile and embedded applications. It introduced depth-wise separable
convolutions and point-wise convolutions to decrease the number of features.
This network helps the architecture to remain lightweight and is easier to
train with the amount of data we have.

###  2.2. Priority maps and biases

Saliency in images exhibits a phenomenon called center bias, where the gazes
of observers are more oriented toward the center of the image compared to the
edges and corners of the stimulus. That is mainly due to the tendency for
relevant semantic information to be located around central regions of images.
The center bias phenomenon can be modeled using a probability map representing
a 2-dimensional Gaussian function as shown below:

(1)
$f(x,y)=\frac{1}{2\pi\sigma_{x}\sigma_{y}}\exp^{-((\frac{x-\mu_{x}}{2\sigma_{x}})+(\frac{y-\mu_{y}}{2\sigma_{y}}))}$

where $x,y$ represent the coordinates of a point on the map.
$(\mu_{x},\mu_{y})$ and $(\sigma_{x},\sigma_{y})$ are the corresponding means
and standard deviations of the distribution, respectively.

In this work, in order to generalize this phenomenon. Our model learns 16 maps
following the same gaussian laws with different means ($\mu$) and standard
deviations ($\sigma$). The integration of this module allows the feature
extractor to only focus on modeling the stimuli-specific features without
handling biases related to the saliency task.

###  2.3. Decoder

The decoder is composed mainly of a merging module, a channel-wise selection
network, and the Soft-ArgMax function. Our merging module network consists of
8 ($3$x$3$) convolutional layers, each activated by a ReLU function. The
succession of these layers ingests the feature maps extracted by the feature
extractor network concatenated with the learnable prior maps previously
described. The merging network combines stimuli-specific features extracted by
MobileNet and prior maps representing the learned biases models, while
gradually reducing the number of feature maps to 20.

An activation probabilities vector from the channel selection network weights
the resulting feature maps. More details about the selection network are given
in section 2.4. The latter allows predicting fixation point coordinates from
the model using the Soft-ArgMax Function (SAM) (Cornia et al., 2018) in a way
similar to (Kerkouri et al., 2021). The SAM function is described as follows:

(2) $SAM(x)=\sum_{i=0}^{W}\sum_{j=0}^{H}\frac{e^{\beta
x_{i,j}}}{\sum_{i^{\prime}=0}^{W}\sum_{j^{\prime}=0}^{H}e^{\beta
x_{i^{\prime},j^{\prime}}}}({\frac{i}{W},\frac{j}{H}})^{T}$

where ${i,j,i^{\prime},j^{\prime}}$ iterate over pixel coordinates. ${H,W}$
represent the height and width of the feature map, respectively. $x$ is the
input feature map and $\beta$ is a parameter adjusting the distribution of the
softmax output.

###  2.4. Channel Selection

The primary function of the channel selection layer is to create a binary mask
vector that chooses the suitable fixation points to be part of the predicted
scanpath. This means it selects the points suitable to be part of the scanpath
and discards other predicted points using a binary mask vector. This selection
mechanism allows generating variable-length scanpaths.

The configuration takes the feature maps originating from the merging module
and passes them to Global Max Pooling and Average Pooling layers in parallel
after flattening. The obtained vectors are used as input to a shallow 3-layer
Multi-Layer Perceptron (MLP) network ( i.e GMPN, GVPN). The resulting vectors
are concatenated, then ingested by a $3^{rd}$ MLP (i.e MN) network activated
by a Sigmoid in the last layer to provide a vector with values in the
continuous range $[0,1]$. This vector is then binarised according to its mean
value as follows:

(3) $mask_{b}=v>mean(v)$

with

(4) $v=MLP_{3}(cat(MLP_{1}(GMaxPool(flat(d))),\\\
MLP_{2}(GAvgPool(flat(d)))))$

where $mask_{b}$ is the resulting binary mask vector.
$MLP_{1},MLP_{2},MLP_{3}$ are the shallow MLP networks, $cat$ is a vector
concatenation, $GMaxPool$ is the global max pooling operation, $GAvgPool$ is
the global average pooling operation, $flat$ is a flattening function and $d$
is the output feature maps obtained from the Merging layers.

Due to the non-differentiability of the binarization function, we cut the
network’s computational graph and multiply the vector $v$ by $d$ to allow the
back-propagation. The Sigmoid function activates the last layer of the
$MLP_{3}$ to avoid an exploding gradient problem and maintain operations in
the same range as the network’s output.

###  2.5. Training and Loss

The model was trained on 9000 images and validated on 1000 from the Salicon
dataset (Jiang et al., 2015). We trained the model using the following loss
function :

(5) $L(y,\hat{y})=BCE(y,\hat{y})+0.001*\sqrt{len(y)^{2}-len(\hat{y})^{2}}$

where $y$ and $\hat{y}$ are the predicted and ground truth scanpath vectors,
respectively. $BCE$ is the Binary Cross-Entropy function and $len$ represents
the length of the scanpath.

The loss function is designed to take into account the difference in length
between scanpaths in addition to the spatial distribution of the fixation
points. We used the Adam optimizer with a learning rate of $5x10^{-5}$ for 70
epochs and initialized all the weights using the Xavier initializer.

###  2.6. Domain Adaptation

The training of our model on natural scene images yields the results shown in
Sec. 3. However, preliminary test results of the same model on images
belonging to other source domains (i.e. paintings, etc.) did not perform as
well. Therefore, to fully utilize the power of our models, we need to adapt
them to a new source domain. The model must be able to perceive the two
domains (i.e. the original natural scene and painting) as being part of the
same data distribution $\mathcal{D}$ for that purpose we use the adaptation
method proposed in (Ganin et al., 2016). That is achieved by minimizing the
perceived distance between the two distributions of natural images
$\mathcal{D}_{n}$ and paintings $\mathcal{D}_{p}$. This means that our feature
extractor should maximize the ability to identify the mutual useful
information between the two domains and discard the ineffective features for
the prediction task.

As we train our models on images from both distributions, we add a small
branch to the network which classifies the images as being from
$\mathcal{D}_{n}$ or $\mathcal{D}_{p}$ in pseudo-labeled manner. We add a
gradient reversal layer (GRL) at the start of this branch, which reverses the
sign of the gradient flow during back-propagation. Eq. 6 defines forward
propagation while Eq. 7 concerns back-propagation. Both equations are given
below:

(6) $GRL(x)=x$ (7) $GRL(\frac{\partial L}{\partial x})=-\frac{\partial
L}{\partial x}$

where $x$ is the input of the layer, and $\frac{\partial L_{d}}{\partial x}$
represents the gradient of the domain loss $L_{d}$ when back propagating
through the network.

While the classifier network learns to discriminate the images from the 2
domains, the reversal of the gradient sign pushes the feature extractor to
find a unified representational space between the domains thus minimizing the
distance between their distributions. This forces the feature extractor to
disregard the domain-specific features and noises and emphasize the mutual
characteristics of the two domains. This can be modeled as the union of the 2
distributions minus the noise distribution of each of the domains:

(8)
$\mathcal{D}=\mathcal{D}_{p}+\mathcal{D}_{n}-(\mathcal{N}_{p}+\mathcal{N}_{n})$

where $\mathcal{D}_{n},\mathcal{D}_{p}$ and $\mathcal{D}$ are defined as
before and $\mathcal{N}_{p}$ and $\mathcal{N}_{n}$ are the specific noise
distributions of the source domain.

For unsupervised domain adaptation training, we used a mix of 2000 Salicon
(Jiang et al., 2015) images and 2000 unlabeled paintings images scraped from
the internet.

##  3\. Experimental Results

Model | Shape | Direction | Length | Position | MM Score | NSS | Congruency
---|---|---|---|---|---|---|---
PathGAN (Assens et al., 2018) | 0.9608 | 0.5698 | 0.9530 | 0.8172 | 0.8252 | -0.2904 | 0.0825
Le Meur (Le Meur and Liu, 2015) | 0.9505 | 0.6231 | 0.9488 | 0.8605 | 0.8457 | 0.8780 | 0.4784
G-Eymol (Zanca et al., 2019) | 0.9338 | 0.6271 | 0.9521 | 0.8967 | 0,8524 | 0.8727 | 0.3449
SALYPATH (Kerkouri et al., 2021) | 0.9659 | 0.6275 | 0.9521 | 0.8965 | 0,8605 | 0.3472 | 0.4572
Our Model | 0.9702 | 0.6173 | 0.9587 | 0.8968 | 0,8607 | 1.0140 | 0.5170
Table 1.  Results of scanpath prediction on Salicon dataset

Model | Shape | Direction | Length | Position | MM Score | NSS | Congruency
---|---|---|---|---|---|---|---
PathGan (Assens et al., 2018) | 0.9237 | 0.5630 | 0.8929 | 0.8124 | 0.7561 | -0.2750 | 0.0209
DCSM (VGG) (Bao and Chen, 2020) | 0.8720 | 0.6420 | 0.8730 | 0.8160 | 0,8007 | - | -
DCSM (ResNet) (Bao and Chen, 2020) | 0.8780 | 0.5890 | 0.8580 | 0.8220 | 0,7868 | - | -
Le Meur (Le Meur and Liu, 2015) | 0.9241 | 0.6378 | 0.9171 | 0.7749 | 0,8135 | 0.8508 | 0.1974
G-Eymol (Zanca et al., 2019) | 0.8885 | 0.5954 | 0.8580 | 0.7800 | 0,7805 | 0.8700 | 0.1105
SALYPATH (Kerkouri et al., 2021) | 0.9363 | 0.6507 | 0.9046 | 0.7983 | 0,8225 | 0.1595 | 0.0916
Our model | 0.9392 | 0.6152 | 0.9100 | 0.8537 | 0.82952 | 0.8888 | 0.2114

Table 2. Results of scanpath prediction on MIT1003.

###  3.1. Datasets

In this study, three datasets have been used to evaluate our method: Salicon
(Jiang et al., 2015), MIT1003 (Judd et al., 2009), and Le Meur paintings (Le
Meur et al., 2020).

Salicon is a large-scale dataset extracted from MS-COCO (Lin et al., 2014)
dataset. Each image is accompanied by the corresponding saliency map and
scanpaths. Here, we used $9000$ images for the training set and $1000$ images
for the validation set. We also used $5000$ images for testing and comparison
with other models. In our work, we tested approximately 250000 scanpaths on
this dataset. The great number of scanpaths ensures the soundness of the
empirical results.

MIT1003 is one of the most well-known datasets for static image saliency, and
it used usually together with the MIT300 dataset (Judd et al., 2012). MIT1003
is composed of 1003 natural scene images depicting objects, individuals, and
scenes in multiple conditions. Each image is provided together with its
saliency map and scanpaths from 15 observers, this constitutes 15045 scanpaths
available for testing. We use the whole dataset essentially to evaluate the
generalization ability of our model (i.e. cross-dataset). Meaning we did not
fine-tune our model on it.

Le Meur paintings (Le Meur et al., 2020) is a saliency dataset for paintings,
which is composed of 150 images for paintings with their saliency maps. This
dataset is employed to evaluate the ability of our model to generalize the
knowledge acquired from natural scene images to the paintings domain and test
the effects of our domain adaptation on this dataset. The dataset does not
provide observers scanpaths, but it is the only visual attention paintings
dataset publicly available.

###  3.2.  Experimental protocol and metrics

We used the same set of images and partitions described in section 3.1 to test
our model and the compared methods. The results of each model were compared to
the ground truth using three metrics: $MultiMatch$ (Dewhurst et al., 2012),
$NSS$ (Peters et al., 2005) and $Congruency$ (Le Meur et al., 2011).

Multimatch compares two vectors representing the scanpaths through 5 different
characteristics (i.e. Shape, Direction, Length, Position and Duration).

NSS (Peters et al., 2005) represents the average normalized saliency value of
a fixation point from the scanpath on the saliency map, while Congruency (Le
Meur et al., 2011) calculates the average ratio of the number of fixations
that fall within a salient region and the global number of fixations.

###  3.3. Results

Model | NSS | Congruency
---|---|---
Our model (without DA) | 1.3620 | 0.4024
Our model (with DA) | 1.5093 | 0.4244

Table 3. Results of scanpath prediction on Le Meur paintings before and after
domain adaptation. Figure 2. Qualitative Results on Salicon Validation. Figure
3. Predictions on Le Meur dataset Before and After Domain Adaptation.

The testing results on Salicon are shown in Table 1. The highest scores are
highlighted in bold. As can be seen, our model achieves the best results on
the $Shape$, $Length$ and $Position$ components while demonstrating acceptable
results for the $Direction$ component. An slightly surpasses the runner-up on
$MM_{S}core$ which represents the mean score of the 4 components. Overall, our
model’s performances on $MultiMatch$ metrics are comparable to the state-of-
the-art models; it even slightly surpasses SALYPATH in some criteria. As far
as it concerns the $NSS$ metrics’ scores, our model achieves the best results
exceeding Le Meur and G-Eymol, which use saliency maps for their scanpath
generation. Yet, the model is surpassed by Le Meur and SALYPATH on the
$congruency$ metric. That can be explained by the fact that Le Meur heavily
relies on saliency maps for scanpath generation while SALYPATH learned
saliency intermediate features from the saliency prediction branch.

Table 2 shows the results obtained on MIT1003. As it can be seen, our model
scores the highest on the $Shape$ and $Position$ criteria while its
performances on the $Direction$ and $Length$ characteristics score slightly
lower albeit still competitive to those of the state-of-the-art methods.
Nonetheless our model surpasses the other models mean score $MMScore$. Our
model $NSS$ and $Congruency$ scores decrease a little. Still they scored
higher than the comparison models. In Figure 2, a qualitative depiction of the
results predicted by our model is given. The figure includes the predicted
scanpaths along with ground truth of scanpaths. These results show that our
model can efficiently and effectively predict plausible scanpaths from images.

###  3.4. Domain Adaptation Results

We evaluated the performance of the domain adaptation of model using a
painting dataset by Le Meur. Unfortunately, the dataset provides only stimuli
images and saliency maps with no scanpaths. That prevented us from using the
$MultiMatch$ metric in our evaluation. Thus we only used the hybrid metrics
$NSS$ and $Congruency$. Table 3 aims at highlighting the improvements of our
model after introducing domain adaption. $NSS$ metrics’ scores increased
considerably and a improvement can even be noticed on the $Congruency$ results
which shows the effectiveness of the domain adaptation in making a closer
distance between the two different datasets. Fig. 3 provides a visual and
qualitative demonstration of the domain adaptation to the scanpath generation
on the new domain of paintings. We can observe that the scanpaths before
domain adaptation have very small saccadic lengths and concentrate the
fixations on a very small region which is sometimes not salient, this can be
attributed to the effect of painting styles. On the other hand the scanpaths
obtained after the domain adaptation span wider regions and is more correlated
to saliency maps depicted in the last row.

##  4\. Conclusion

In this paper, we presented a new deep learning architecture to predict
scanpath for 2D images, we later adapt the approach to painting images which
represent a fundamental segment of cultural heritage. Our model consists of a
feature extractor reinforced task oriented learnable prior maps. The network
also uses a channel-wise selection module that ensures the stochastic nature
of the scanpath length between images. The proposed model was evaluated
against state-of-the-art ones and produced the best results. Using an
adversarial unsupervised domain adaptation algorithm, the model has later
adapted from the natural scenes domain to the paintings domain, producing very
relevant results. The Soft-ArgMax function creates limitation as it is not
able to predict the duration of fixations, in future works we will add a time
prediction module as well as a incorporate a probabilistic generative approach
in order to make the prediction process of stochastic nature. Nonetheless
incorporating the domain adaptation approach for scanpath prediction, opens
doors for to further expand the model to other image domains related to
heritage understanding like monuments architectures, statues, and artificial
3D images applications like virtual museum visits.

###### Acknowledgements.

Funded by the TIC-ART project, Regional fund (Region Centre-Val de Loire)

## References

  * (1)
  * Abouelaziz et al. (2019) Ilyass Abouelaziz, Aladine Chetouani, Mohammed El Hassouni, Lj Latecki, and Hocine Cherifi. 2019\.  3D visual saliency and convolutional neural network for blind mesh quality assessment.  _Neural Computing and Applications_ (2019).  https://doi.org/10.1007/s00521-019-04521-1(
  * Abouelaziz et al. (2020) Ilyass Abouelaziz, Aladine Chetouani, Mohammed El Hassouni, Longin Jan Latecki, and Hocine Cherifi. 2020.  No-reference mesh visual quality assessment via ensemble of convolutional neural networks and compact multi-linear pooling.  _Pattern Recognition_ 100 (2020), 107174.  https://doi.org/10.1016/j.patcog.2019.107174
  * Assens et al. (2018) Marc Assens, Xavier Giro-i Nieto, Kevin McGuinness, and Noel E O’Connor. 2018.  PathGAN: visual scanpath prediction with generative adversarial networks. In _Proceedings of the European Conference on Computer Vision (ECCV)_. 0–0.
  * Bao and Chen (2020) Wentao Bao and Zhenzhong Chen. 2020.  Human Scanpath Prediction based on Deep Convolutional Saccadic Model.  _Neurocomputing_ (2020).
  * Bruno et al. (2019) Alessandro Bruno, Francesco Gugliuzza, Edoardo Ardizzone, Calogero Carlo Giunta, and Roberto Pirrone. 2019.  Image content enhancement through salient regions segmentation for people with color vision deficiencies.  _i-Perception_ 10, 3 (2019), 2041669519841073\.
  * Chetouani (2018a) Aladine Chetouani. 2018a.  A Blind Image Quality Metric using a Selection of Relevant Patches based on Convolutional Neural Network. In _EUSIPCO_. 1452–1456.  https://doi.org/10.23919/EUSIPCO.2018.8553127
  * Chetouani (2018b) Aladine Chetouani. 2018b.  Convolutional neural network and saliency selection for blind image quality assessment. In _2018 25th IEEE International Conference on Image Processing (ICIP)_. IEEE, 2835–2839.
  * Chetouani and Li (2020) Aladine Chetouani and Leida Li. 2020.  On the use of a scanpath predictor and convolutional neural network for blind image quality assessment.  _Signal Processing: Image Communication_ 89 (2020), 115963.  https://doi.org/10.1016/j.image.2020.115963
  * Cornia et al. (2017) Marcella Cornia, Lorenzo Baraldi, Giuseppe Serra, and Rita Cucchiara. 2017.  Visual saliency for image captioning in new multimedia services. In _2017 IEEE International Conference on Multimedia & Expo Workshops (ICMEW)_. IEEE, 309–314.
  * Cornia et al. (2018) Marcella Cornia, Lorenzo Baraldi, Giuseppe Serra, and Rita Cucchiara. 2018.  Predicting human eye fixations via an lstm-based saliency attentive model.  _IEEE Transactions on Image Processing_ 27, 10 (2018), 5142–5154.
  * Dewhurst et al. (2012) Richard Dewhurst, Marcus Nyström, Halszka Jarodzka, Tom Foulsham, Roger Johansson, and Kenneth Holmqvist. 2012.  It depends on how you look at it: Scanpath comparison in multiple dimensions with MultiMatch, a vector-based approach.  _Behavior research methods_ 44, 4 (2012), 1079–1100.
  * et al (2019) M Hamidi et al. 2019\.  Blind Robust 3D Mesh Watermarking Based on Mesh Saliency and Wavelet Transform for Copyright Protection.  _Information_ 10, 2 (2019).  https://doi.org/10.3390/info10020067
  * Ganin et al. (2016) Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario March, and Victor Lempitsky. 2016.  Domain-Adversarial Training of Neural Networks.  _Journal of Machine Learning Research_ 17, 59 (2016), 1–35.  http://jmlr.org/papers/v17/15-239.html
  * Goldberg and Helfman (2010) Joseph H. Goldberg and Jonathan I. Helfman. 2010.  Visual Scanpath Representation. In _Proceedings of the 2010 Symposium on Eye-Tracking Research & Applications_ (Austin, Texas) _(ETRA ’10)_. Association for Computing Machinery, New York, NY, USA, 203–210.  https://doi.org/10.1145/1743666.1743717
  * Hartmann (2020) Rudi Hartmann. 2020\.  Heritage and Economy.  In _International Encyclopedia of Human Geography (Second Edition)_ (second edition ed.), Audrey Kobayashi (Ed.). Elsevier, Oxford, 369–372.  https://doi.org/10.1016/B978-0-08-102295-5.10469-X
  * Itti and Koch (2001) Laurent Itti and Christof Koch. 2001.  Computational modelling of visual attention.  _Nature reviews neuroscience_ 2, 3 (2001), 194–203.
  * Jiang et al. (2015) Ming Jiang, Shengsheng Huang, Juanyong Duan, and Qi Zhao. 2015\.  SALICON: Saliency in Context.. In _CVPR_. IEEE Computer Society, 1072–1080.  http://dblp.uni-trier.de/db/conf/cvpr/cvpr2015.html#JiangHDZ15
  * Judd et al. (2012) Tilke Judd, Frédo Durand, and Antonio Torralba. 2012\.  A Benchmark of Computational Models of Saliency to Predict Human Fixations. In _MIT Technical Report_.
  * Judd et al. (2009) Tilke Judd, Krista Ehinger, Frédo Durand, and Antonio Torralba. 2009.  Learning to Predict Where Humans Look. In _IEEE International Conference on Computer Vision (ICCV)_.
  * Kerkouri and Chetouani (2021) Mohamed Amine Kerkouri and Aladine Chetouani. 2021.  A Simple and efficient deep Scanpath Prediction.  arXiv:2112.04610 [cs.CV]
  * Kerkouri et al. (2021) Mohamed A. Kerkouri, Marouane Tliba, Aladine Chetouani, and Rachid Harba. 2021.  Salypath: A Deep-Based Architecture For Visual Attention Prediction. In _2021 IEEE International Conference on Image Processing (ICIP)_. 1464–1468.  https://doi.org/10.1109/ICIP42928.2021.9506295
  * Koch and Ullman (1987) Christof Koch and Shimon Ullman. 1987.  Shifts in selective visual attention: towards the underlying neural circuitry.  In _Matters of intelligence_. Springer, 115–141.
  * Le Meur et al. (2011) Olivier Le Meur, Thierry Baccino, and Aline Roumy. 2011\.  Prediction of the inter-observer visual congruency (IOVC) and application to image ranking. In _Proceedings of the 19th ACM international conference on Multimedia_. 373–382.
  * Le Meur et al. (2020) Olivier Le Meur, Tugdual Le Pen, and Rémi Cozot. 2020\.  Can we accurately predict where we look at paintings?  _PLOS ONE_ 15 (10 2020), 1–20.  https://doi.org/10.1371/journal.pone.0239980
  * Le Meur and Liu (2015) Olivier Le Meur and Zhi Liu. 2015.  Saccadic model of eye movements for free-viewing condition.  _Vision research_ 116 (2015), 152–164.
  * Lin et al. (2014) Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. 2014.  Microsoft coco: Common objects in context. In _European conference on computer vision_. Springer, 740–755.
  * Mahmoudi et al. (2020a) M.Amine Mahmoudi, Aladine Chetouani, Fatma Boufera, and Hedi Tabia. 2020a.  Kernelized Dense Layers For Facial Expression Recognition. In _ICIP)_. 2226–2230.  https://doi.org/10.1109/ICIP40778.2020.9190694
  * Mahmoudi et al. (2020b) M. Amine Mahmoudi, Aladine Chetouani, Fatma Boufera, and Hedi Tabia. 2020b.  Learnable pooling weights for facial expression recognition.  _Pattern Recognition Letters_ 138 (2020), 644–650.  https://doi.org/10.1016/j.patrec.2020.09.001
  * Patel et al. (2021) Yash Patel, Srikar Appalaraju, and R Manmatha. 2021\.  Saliency Driven Perceptual Image Compression. In _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_. 227–236.
  * Peters et al. (2005) Robert J Peters, Asha Iyer, Laurent Itti, and Christof Koch. 2005.  Components of bottom-up gaze allocation in natural images.  _Vision research_ 45, 18 (2005), 2397–2416.
  * Raichle (2010) Marcus E Raichle. 2010\.  The brain’s dark energy.  _Scientific American_ 302, 3 (2010), 44–49.
  * Sandler et al. (2018) Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen. 2018.  Mobilenetv2: Inverted residuals and linear bottlenecks. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 4510–4520.
  * Tliba et al. (2022a) Marouane Tliba et al. 2022a.  SATSal: A Multi-level Self-Attention Based Architecture for Visual Saliency Prediction.  _IEEE Access_ (2022), 1–1.  https://doi.org/10.1109/ACCESS.2022.3152189
  * Tliba et al. (2022b) Marouane Tliba, Mohamed Amine Kerkouri, Aladine Chetouani, and Alessandro Bruno. 2022b.  Self Supervised Scanpath Prediction Framework for Painting Images. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_. 1539–1548.
  * Treisman and Gelade (1980) Anne M Treisman and Garry Gelade. 1980.  A feature-integration theory of attention.  _Cognitive psychology_ 12, 1 (1980), 97–136.
  * Wang et al. (2020) Haoxiang Wang, Zhihui Li, Yang Li, Brij B Gupta, and Chang Choi. 2020.  Visual saliency guided complex image retrieval.  _Pattern Recognition Letters_ 130 (2020), 64–72.
  * Yue and Jun (2021) Yin Yue and Zhou Jun. 2021\.  Rock painting restoration method based on saliency detection and TV model. In _2021 International Conference on Computer Engineering and Application (ICCEA)_. 122–125.  https://doi.org/10.1109/ICCEA53728.2021.00032
  * Zanca et al. (2019) Dario Zanca, Stefano Melacci, and Marco Gori. 2019\.  Gravitational laws of focus of attention.  _IEEE transactions on pattern analysis and machine intelligence_ (2019).


# Tunable drag drop via flow-induced snap-through in origami

Rishabh Nain  Tom Marzin  Laboratoire d’Hydrodynamique, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Chemical and Biological Engineering, Princeton University, NJ 08544, United
States  Sophie Ramananarivo  Laboratoire d’Hydrodynamique, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

###### Abstract

We leverage the snap-through response of a bistable origami mechanism to
induce a discontinuous evolution of drag with flow speed. The transition
between equilibrium states is passively actuated by airflow, and we
demonstrate that large shape reconfiguration over a small increment of flow
velocity leads to a pronounced and sudden drop in drag. Moreover, we show that
systematically varying the geometrical and mechanical properties of the
origami unit enables the tuning of this drag discontinuity and the critical
speed and loading at which it occurs. Experimental results are supported by a
theoretical aero-elastic model, which further guides inverse design to
identify the combination of structural origami parameters for targeted drag
collapse. This approach sheds light on harnessing origami-inspired mechanisms
for efficient passive drag control in a fluid environment, applicable for load
alleviation or situations requiring swift transitions in aerodynamic
performances.

The drag force acting on an object typically increases with the speed of the
surrounding fluid flow, as the pressure and shear stress exerted on its
surface grow. For flexible bodies, this increase can be partially offset by
the deformations induced by fluid forces, as it produces a more streamlined
shape and reduced frontal area [1, 2, 3, 4, 5, 6]. This reconfiguration
reduces the load compared to that of a non-deformable body; nonetheless, drag
still monotonically increases with speed, although at a lower rate. However,
there are a few instances where drag has been observed to drop sharply with
velocity. One such phenomenon occurs due to a sudden change in flow regime,
referred to as drag crisis, which happens when the boundary layer transitions
from laminar to turbulent, moving its separation downstream and narrowing the
wake [7, 8]. Another scenario involves abrupt and significant shape changes in
the body, which can notably occur due to partial breakage, as observed in
wind-induced tree pruning [9]. While this process effectively provokes a sharp
decrease in the drag force and preserves the integrity of the rest of the
structure, it is not reversible.

An attractive alternative is snap-through events in flexible structures,
during which the elastic system undergoes a rapid transition from one
equilibrium state to another. This process is swift, repeatable and produces
significant shape changes, making it increasingly popular in engineering for
obtaining switchable optical properties, reusable energy trapping materials,
microfluidic pumps, or in soft robotics [10, 11, 12, 13]. Snapping
instabilities occur through the application of an external input like a point
load, electrostatic forces, or thermal effects, and can also be triggered by
fluid flows. For example, a post-buckled sheet snaps-through and snaps-back in
response to air flowing at different speeds, thus opening and closing an air
inlet for flow regulation applications [14, 15]. The snapping of a slender
membrane or arc has also been used to modify a channel resistance for passive
control of viscous flow [16, 17], or for energy harvesting from wind-induced
snapping oscillations [18, 19]. The critical flow speed to initiate the
instability can be adjusted through the non-linear structural mechanics of the
post-buckled member. Note that the system passively responds to flows, which
makes it inherently autonomous and leverages the need for additional complex
actuation. A few studies focused on the effect of such fast-changing
geometries on aerodynamic loading. In the context of wind turbines, bistable
winglets or airfoil trailing-edges autonomously transition to another stable
state that alleviates excessive fluid load [20, 21]. While reported drag
reduction is limited to about 15% [21], primarily due to a focus on lift
applications, it shows potential for more substantial drag drops with
appropriate designs.

A promising framework for manufacturing shape-morphing structures with non-
linear mechanical properties is the technique of origami, inspired by the
Japanese art of paper folding. A surface is equipped with a set of folds,
which act as a linkage mechanism, producing large three-dimensional shape
changes along selected kinematic pathways that are prescribed by the folding
pattern. Origami has also proven to be a powerful method to tailor targeted
mechanical properties, and in particular multistability [22]. Some specimens
can thus settle in more than one stable state in the absence of external load
and produce a rapid snap-through response as they pass from one equilibrium to
the other [23, 24, 25, 22]. It offers predictable and reversible shape
transformation, with the possibility to tune the shape and number of
designated stable configurations through the fold arrangement, crease
torsional stiffness, and assigned rest angle [26]. Those transformable
structures with programmable snap-through mechanisms have the potential to
produce a controlled drop in drag at targeted flow conditions.

As a starting point, we demonstrate the potential for drag collapse with the
waterbomb base, a single vertex bistable origami mechanism that is one of the
fundamental origami units [23, 27]. This umbrella-like cell has two stable
states and features a rapid snap-through response as it transitions from one
to the other. Origami structures have already raised interest for drag
control, due to their large morphing capabilities [28, 29, 30, 31]. In
particular, a previous study on the waterbomb base showed that its extreme
reconfiguration tends to reduce drag to the point that it no longer increases
with flow speed in the regime of large deformation [30]. Here, we make use of
its bistable nature, which has not been utilized before for drag applications.
We show that it produces a discontinuous evolution of drag with flow speed,
marked by a sudden and significant drop as the structure snaps through. We
systematically vary the geometrical and mechanical properties of the unit and
show that it allows us to tune the drag discontinuity and the critical speed
and loading at which it occurs. Experiments are captured by a simplified
theoretical model that is further used for inverse design, which is finding
the structural origami parameters that produce a targeted drag collapse.

##  1 Results

Drag collapse through snapping

Figure 1: Drag collapse through flow-induced snap-through of a bistable
origami unit. (a) Waterbomb base with radius $R$, composed of eight
alternating valley and mountain folds with respective angles $\psi_{v}$ and
$\psi_{m}$, and with unimodal kinematics (for symmetric folding) described by
the angle $\theta$. In the absence of loading, the unit sits in one of two
stable equilibria with angle $\theta_{0}$ and $\theta_{0}^{*}$. Evolution of
(b) the drag force $F_{d}$ on the bistable origami cell when exposed to a
uniform airflow with increasing velocity $U$, and of (c) its dimensionless
frontal area $S/\pi R^{2}$ projected in the plane perpendicular to the flow
(see schematics of the experimental set-up in the inset of b). The drag
collapse observed in the snapping unit (with $R=5$ cm, $\kappa=19\times
10^{-3}$ N and $\theta_{0}=41^{\circ}$, black solid curves) is contrasted with
the quadratic drag experienced by a rigid unit with the same geometry (same
$R$ and $\theta_{0}$, grey dashed curves).

The waterbomb base consists of a disk of radius $R$, featuring alternating
mountain and valley folds distributed evenly around a central vertex. The
present study focuses on the traditional version of the origami waterbomb,
comprising a total of eight folds, and which is rigid-foldable. Namely, facets
are considered rigid, with all movement occurring at the joints, thus
simplifying the kinematic analysis. Additionally, we assume symmetry of
motion, based on experimental observation, meaning that all the mountain folds
deflect by the same amount, and the same holds for the valley folds. As a
result, the mechanism reduces to one degree-of-freedom system [23, 32], where
the shape of the entire structure can be determined from a single variable,
chosen as the angle $\theta$ between a valley fold and the central axis
passing through vertex (see Fig.1a).

Importantly, this unimodal origami mechanism has two stable equilibrium states
and can transition between them through snapping. When at rest, the folded
unit sits in one of these stable positions with its vertex pointing in one
direction (see Fig.1a-top). Upon applying an external force, the structure
flattens out towards a plane that represents an unstable equilibrium. When
displaced beyond this plane, it snaps through and moves to a second stable
position with the vertex pointing in the opposite direction (see
Fig.1a-bottom). These two stable positions are characterized by angles
$\theta_{0}$ and $\theta_{0}^{*}$, and although they may look similar in
Fig.1a, they are not mirror images of each other. The relationship between the
angles of mountain and valley folds, $\Psi_{m}$ and $\Psi_{v}$ (measured in
terms of deviation of straightness, see Fig.1a), and $\theta$ also differs on
both sides of the unstable equilibrium [23].

The origami bistable mechanism presents the advantage of being readily
tunable. The stiffness of the folds notably influences the forces required to
transition between stable states. One additional tuning parameter is the rest
configuration, characterized by the angle $\theta_{0}$ in the first stable
state. A compact unit (with a smaller $\theta_{0}$) requires more elastic
energy to pass through the unstable planar state, as it has to undergo a
greater deflection. The second equilibrium state $\theta_{0}^{*}$ also changes
correspondingly with $\theta_{0}$, modulating the amplitude of the shape
change upon snapping (the relationship between $\theta_{0}$ and
$\theta_{0}^{*}$ is illustrated in Supplementary Fig.5).

Specimens are fabricated using a technique developed in a previous study [30]
and inspired by the literature [33, 34], which is based on layer
superposition. The rigid facets are laser cut from a 350$\mu$m thick sheet of
Mylar, which is then attached to a thinner sheet using double-sided tape,
forming the creases. These folds act as elastic hinges of specific stiffness
and rest angles. Rest angles are prescribed by utilizing the thermo-formable
property of PET (see Methods). As shown in a previous study, the folds then
operate as torsional springs, exhibiting elastic behavior around their
assigned rest angle [30]. The stiffness can be adjusted by using different
thicknesses of the thin sheet, ranging from $40-100\mu$m. To measure the
stiffness experimentally, a single fold is fabricated using the same process
and is tested by tracking the opening angle as an increasing torque is applied
to the hinge (see Methods and Supplementary Fig.1). The value of the stiffness
is then obtained from the slope of the torque-angle curve. The fabrication
process shows a standard deviation of 8$\%$ for $\kappa$.

We conducted tests on ten origami units, sequentially varying the radius
$R=4-7$ cm, the stiffness of the folds $\kappa=8-92\times 10^{-3}$ N, and the
rest configuration characterized by the angle
$\theta_{0}=32^{\circ}-66^{\circ}$. These units are subjected to a steady air
flow generated by an open jet wind tunnel with a square cross-section of width
$40$ cm and flow speed varied in $1-16$ m/s. An elbow arm is used upstream, to
hold the unit at its vertex in the center of the cross-section (see inset of
Fig.1b). This arm is connected to a six-component force sensor, measuring the
drag component of the fluid force in the direction of the flow. For each
measurement, we average the drag over 30 seconds and subtract the drag on the
support system alone. The flow-induced shape changes are characterized through
the projected area $S$ of the unit in the plane perpendicular to the flow,
captured by a camera placed downstream and extracted through Matlab custom
image analysis. The observed kinematics align with the unimodal mechanism
described earlier, which is rigidly foldable with motion restricted to the
joint and exhibits rotational symmetry of motion. The entire geometry of the
unit can be reconstructed from $S$, which will be further used as the
observable to measure shape reconfiguration (the relationship between $S$ and
$\theta$ is illustrated in Supplementary Fig.4).

Fig.1b-c displays the results of a typical experiment. The unit is initially
in its first stable position $\theta_{0}$, with the concave side facing the
flow. As the flow speed gradually increases, the cell symmetrically opens up.
The expansion of the frontal area, reported in Fig.1c, leads to a faster
increase of drag with flow speed (solid black curve in Fig.1b) compared to the
classical $U^{2}$ law observed for a rigid cell frozen in its initial
configuration (dashed grey curve in Fig.1b). When the cell reaches the
unstable flat state, it snaps to a more compact and streamlined configuration.
Note that slight cell vibrations make it challenging to precisely attain
$S/\pi R^{2}=1$. The sudden reconfiguration results in a discontinuity in the
drag force on the unit, and the shape change is significant enough to produce
a drag drop by up to $69\%$ among the specimen tested in this study. As $U$ is
further increased, the unit closes up, reducing the frontal area exposed to
flow. This closing regime has been extensively studied in a previous study
[30]; it features an increase of drag that is slower than the quadratic
$U^{2}$ law for rigid objects, owing to the reduction in frontal area combined
with shape streamlining. Variability in the experiments mainly results from
the origami fabrication process, yielding a relative dispersion of about $7\%$
for the drag data and $5\%$ for the critical snapping velocity (see
Supplementary Fig.3).

We thus observe a non-continuous evolution of drag with flow speed with a
marked collapse. This behavior can be described by selected observable
parameters, namely the critical flow velocity of the snapping $U_{c}$, the
maximum drag force reached just before snapping $F_{d,max}$, and the
subsequent drag drop $\Delta F_{d}$. These parameters are prominent features
characterizing the discontinuity but are also relevant for practical
applications, such as designing protection devices or valves. They determine a
maximum load, the amount of load reduction, or the critical speed at which the
structure would operate. In the following section, we investigate how those
quantities are related to the origami geometrical and mechanical properties,
specifically its size $R$, rest configuration $\theta_{0}$, and fold stiffness
$\kappa$.

Aeroelastic Mechanism

Figure 2: Influence of the origami structural parameters on drag. Evolution of
(a, c, and e) the drag force $F_{d}$ and (b, d, and f) the dimensionless
frontal area $S/\pi R^{2}$ as a function of flow velocity $U$, for origami
units with different (a-b) radii $R$, (c-d) fold stiffness $\kappa$, and (e-f)
rest angle $\theta_{0}$. The unit $R=5$ cm, $\kappa=19\times 10^{-3}N$ and
$\theta_{0}=41^{\circ}$ is common to all graphs.

We first vary the radius $R$ of the unit, while keeping the stiffness to be
$\kappa=19\times 10^{-3}$ N and the rest angle to be $\theta_{0}=41^{\circ}$.
Note that fabrication produces a small variability in the rest angle among
specimens, which is visible in Fig.2b in the $6\%$ dispersion in the
dimensionless frontal area in the absence of flow ($U=0$). As shown in
Fig.2a-b, all units exhibit the characteristic behavior described earlier,
that is a non-continuous evolution of drag with flow speed with a sudden drop
associated with the snap-through of the origami unit. However, bigger units
(darker orange curves) initially experience a higher fluid loading, because of
the larger area exposed to the flow. As a result, they reach the flat state
$S/\pi R^{2}=1$ faster and flip at a lower critical speed (see Fig.2b).
Interestingly, in Fig.2a, the maximum drag $F_{d,max}$ experienced by the
units before snapping is the same, and they feature similar jumps amplitude
$\Delta F_{d}$. Changing the radius of the unit thus primarily shifts the
critical velocity $U_{c}$ at which the drag collapse occurs, while preserving
$F_{d,max}$ and $\Delta F_{d}$.

Fig.2c-d shows the results for origami cells with different fold stiffness,
but the same size $R=5$ cm and rest angle $\theta_{0}=41^{\circ}$ (with a
similar degree of variability as earlier). Stiffer units (darker blue curves)
deploy slower with increasing $U$ owing to larger resistance to deformation,
and reach the flat state at larger critical velocities (see Fig.2d). The most
rigid one does not reach this unstable point within the range of flow covered
here. Snapping at larger speeds $U_{c}$ results both in higher drag peak
values $F_{d,max}$ and larger jumps $\Delta F_{d}$ (see Fig.2c). Stiffness
thus affects all three drag features, offering an additional lever for
control.

Finally, we vary the degree of opening of the cell rest state $\theta_{0}$,
while keeping the size and stiffness constant as $R=5$ cm and $\kappa=19\times
10^{-3}$ N respectively. Starting from a more compact rest state (lighter
green curves) delays the flipping towards higher flow velocities, as depicted
in Fig.2f. Intuitively, reaching the flat state that is further away requires
greater work from fluid loading. As was the case for the $\kappa$-series of
units, snapping at greater $U_{c}$ results in larger maximal drag $F_{d,max}$
in Fig.2e. As mentioned earlier and illustrated in the schematics of Fig.2f,
changing $\theta_{0}$ also affects the degree of closure of the second stable
state. Small $\theta_{0}$ (light green) results in a larger collapse in the
frontal area post-snapping compared to units initially more opened (dark
green) in Fig.2f. It correspondingly leads to larger jumps in drag $\Delta
F_{d}$ in Fig.2e. The rest angle thus completes the set of three parameters
that will allow us to adjust the three drag characteristics $U_{c}$,
$F_{d,max}$, and $\Delta F_{d}$.

Figure 3: Aeroelastic mechanism. Evolution of the dimensionless drag
$\overline{F_{d}}=F_{d}/\kappa$ as a function of the Cauchy number $C_{y}=\rho
U^{2}R^{2}/\kappa$ for origami units with (a) different radii $R$ and fold
stiffness $\kappa$ (while maintaining the same rest angle
$\theta_{0}=41^{o}$), and (b) different initial degree of opening $\theta_{0}$
(while keeping the same $R=5$ cm and $\kappa=19\times 10^{-3}$ N).
Experimental data are presented with the same markers and color code as in
Fig.2, and compared to theoretical predictions (solid curves). In panel b,
theoretical curves are displayed in a separate graph for readability, and with
the same color gradient as in experiments indicating different $\theta_{0}$.

The trends observed through variations in parameters reflect the interplay
between external fluid loading and elastic restoring forces of the origami
unit, which governs the unit’s reconfiguration process. This fluid-elastic
competition can be analyzed in terms of Cauchy number $C_{y}$, as introduced
in a previous study [30] and frequently used in the literature [2, 3, 4]. This
dimensionless number examines the relative magnitudes of the work done by
fluid forces, $\rho U^{2}R^{3}$, and the elastic energy of folding, $\kappa
R$:

$C_{y}=\frac{\rho U^{2}R^{2}}{\kappa}$ (1)

In Fig.3a, drag measurements are re-plotted in the dimensionless form
$\overline{F_{d}}=F_{d}/\kappa$ as a function of $C_{y}$, using the fold
stiffness $\kappa$ as a characteristic force scale, in line with [30]. Here,
we only report data for units with different stiffness $\kappa$ and radius
$R$, which share the same initial rest configuration $\theta_{0}$. Data
collapse onto a master curve, with a drag peak and subsequent jump that are
solely set by the fold stiffness, and which occur at a critical Cauchy number
of $C_{y,c}\approx 14$ for all units. The dimensionless quantities $C_{y}$ and
$\overline{F_{d}}$ capture the effect of origami size and stiffness in fluid-
elastic mechanisms. However, it does not account for the more intricate
influence of the rest angle $\theta_{0}$. As shown in the top panel of Fig.3b,
the latter impacts the peak and jump value of $\overline{F_{d}}$, as well as
the critical Cauchy number for snapping.

Next, we develop a simplified theoretical model to further explore and
corroborate the link between the origami unit’s characteristics and drag. This
model will also be used for the inverse problem. It employs the same energy
approach as previous work on the closing dynamics of the waterbomb unit in a
flow [30] and has been extended to account for the snapping. We provide a
brief overview of the main elements of the model here, and additional details
can be found in the Methods and in [30]. The static equilibrium of the unit in
a flow is determined through the minimization of its energy $E_{el}-W$,
comprising both the elastic potential energy associated with folding $E_{el}$
and the work done by fluid forces $W$ as the unit deforms:

$\displaystyle E_{el}-W=4\frac{\kappa
R}{2}[(\psi_{m}-\psi_{m}^{0})^{2}+(\psi_{v}-\psi_{v}^{0})^{2}]$ (2)
$\displaystyle-8\int_{\theta_{0}}^{\theta}\int_{S_{f}}\rho(\mathbf{U}.\mathbf{n}_{\theta^{\prime}})^{2}dS_{f}\mathbf{n}_{\theta^{\prime}}.d\mathbf{X}_{\theta^{\prime}}$

Where $\rho$ is the air density and $S_{f}=\pi R^{2}/8$ is the facet area
(with other quantities introduced subsequently). The first term of Eq.2
corresponds to the elastic potential energy of the four mountain and four
valley folds, modeled as elastic torsional springs with stiffness $\kappa$. It
scales with the square of the deviation of their respective angles $\psi_{m}$
and $\psi_{v}$ to their rest values $\psi_{m}^{0}$ and $\psi_{v}^{0}$ [23].
The second term corresponds to the work done by fluid forces acting on the
unit facets, which is the same for all eight facets owing to the rotational
symmetry of motion. The Reynolds number for our experiments falls within the
range $Re=UR/\nu\approx 10^{3}-10^{5}$, where form drag predominates over
friction drag, and the latter is neglected. As discussed in [30], friction
forces likely become significant in regimes with substantial deformation,
where facets nearly align with the flow. However, this particular regime is
not reached in our present study. Consistent with earlier studies on the
reconfiguration of flexible structures in a flow [3, 5, 6, 30], pressure is
computed using conservation of momentum principles. It is proportional to the
square of the component of flow momentum perpendicular to the facet, denoted
as $\rho(\mathbf{U}\cdot\mathbf{n}_{\theta})^{2}$, where $\mathbf{n}_{\theta}$
represents the normal unit vector. While simplified, this model has
demonstrated effectiveness in capturing fluid-elastic behaviors [3, 5, 6, 30].
It offers a reasonable approximation of more complex potential flow
expressions [35, 30], which would be challenging to implement for three-
dimensional geometries. The pressure force is integrated over the facet area
$S_{f}$ and its work is computed along the facet trajectory, with
$d\mathbf{X}_{\theta}$ the local displacement. By expressing the energy in a
dimensionless form as $\overline{E_{el}-W}=(E_{el}-W)/\kappa R$, the
previously introduced Cauchy number $C_{y}$ emerges as the pertinent fluid-
elastic parameter:

$\displaystyle\overline{E_{el}-W}=2[(\psi_{m}-\psi_{m}^{0})^{2}+(\psi_{v}-\psi_{v}^{0})^{2}]$
(3)
$\displaystyle-8C_{y}\int_{\theta_{0}}^{\theta}\int_{\overline{S}_{f}}\rho(\mathbf{e}_{z}.\mathbf{n}_{\theta^{\prime}})^{2}d\overline{S}_{f}\mathbf{n}_{\theta^{\prime}}.d\overline{\mathbf{X}}_{\theta^{\prime}}$

Where $\bar{.}$ denotes dimensionless quantities, using the radius as
characteristic length, and $\mathbf{e}_{z}$ is the unit vector in the
direction of the flow. In Eq. 3, all kinematic variables of the 1-DOF
mechanism depend on a single parameter $\theta$. However, these relationships
differ on each side of the flat state. The equilibrium angle $\theta$, which
satisfies $d\overline{(E_{el}-W)}/d\theta=0$, is determined numerically while
gradually increasing the Cauchy number, first addressing the opening phase of
the kinematics. As the system approaches the flat state at $\theta=90^{o}$
beyond a specific criterion, the kinematic relations are modified to describe
the closing phase. From the static equilibrium angle $\theta$, we can
determine the drag force acting on the unit as the component of fluid pressure
forces in the direction of the flow (see Methods).

Theoretical drag results are presented in Fig.3a in dimensionless form,
considering a given rest angle $\theta_{0}=41^{o}$ based on experiments (black
solid line). The simplified model successfully captures experimental features,
displaying a peak in drag that remains unaffected by the unit’s size, followed
by a sudden drop at a critical Cauchy number. Although the model tends to
overestimate $F_{d,max}$ and $\Delta F_{d}$, the quantitative agreement is
reasonable, considering that no adjustable parameters were used and the
simplifications made. As depicted in the lower panel of Fig.3b, the model also
effectively captures the impact of varying the initial rest angle $\theta_{0}$
on drag. A more compact unit exhibits a higher drag peak and jump, along with
a higher critical Cauchy number. This model will be valuable for the inverse
design of origami units to achieve the desired drag, as we will discuss in the
next section.

Inverse design approach

Figure 4: Targeting drag collapse features. (a) The inverse design approach
consists of finding the set of origami structural parameters ($R$, $\kappa$,
$\theta_{0}$) that will lead to a collapse of drag with a targeted maximum
before snapping $F_{d,max}$, jump $\Delta F_{d}$, and occurring at a defined
critical flow speed $U_{c}$. Evolution of drag with flow speed for physical
prototypes of optimized origami designs, which are calculated to produce a
selective variation of (b) $U_{c}$, (c) $\Delta F_{d}$, and (d) $F_{d,max}$,
while keeping the two remaining features identical. Respective target values
of ($U_{c}$, $\Delta F_{d}$, $F_{d,max}$) are indicated in the Methods.

We have characterized how the waterbomb unit features influence its drag
behavior. We now pose the question: “Can we identify the specific geometrical
and mechanical origami parameters that would lead to a targeted drag collapse
in a uniform fluid flow?” More specifically, we formulate the inverse problem
where the objective is to find the optimal set of parameters
$(R,\kappa,\theta_{0})$ that will result in a collapse of drag, with a
predefined target peak drag force $F_{d,max}^{t}$ and jump $\Delta F_{d}^{t}$
value at a critical flow speed $U_{c}^{t}$. To address this, we use the model
presented earlier and validated experimentally, which establishes a
relationship between the parameters $(R,\kappa,\theta_{0})$ and the resulting
drag characteristics $(F_{d,max},\Delta F_{d},U_{c})$. The optimization is
implemented numerically in Matlab using the trust-region-dogleg algorithm of
fsolve solver. It solves the system of equations $F_{d}-F_{d,max}^{t}=0$,
$\Delta F_{d}-\Delta F_{d}^{t}=0$ and $U_{c}-U_{c}^{t}=0$ with respect to the
origami unit’s features $(R,\kappa,\theta_{0})$. The algorithm finds a single
solution, satisfying the constraints by a maximum residual of $O(10^{-19})$
within a hundred iterations.

To validate the inverse design approach, we test our prediction experimentally
by fabricating the optimized design and measuring the evolution of its drag
with flow speed. While the specific radius $R$ and rest angle $\theta_{0}$ can
be easily implemented through laser-cutting and an adjustable mount, achieving
an arbitrary fold stiffness $\kappa$ is more challenging. This limitation
arises from the finite number of sheet thicknesses available for manufacturing
the folds. To overcome this constraint, we introduce slits in the folds. These
slits modulate the effective length of the folds and, consequently, their
stiffness, enabling a controlled variation of $\kappa$ by adjusting their
number (see Methods and Supplementary Fig.2).

To demonstrate our ability to control the drag behavior, we create sets of
specimens with one parameter gradually changing while keeping the other two
fixed. We first vary the targeted critical flow speed for snapping, as
$U_{c}^{t}=[4,7,10]$ m/s, while prescribing constant peak drag value
$F_{d,max}^{t}=0.8$ N and jump $\Delta F_{d}^{t}=0.6$ N. The corresponding
optimized origami parameters $(R,\kappa,\theta_{0})$ are provided in the
Methods. Consistent with the differences observed earlier in Fig.3 between
experiments and theory, the experimental realizations exhibit lower values in
Fig.4b, namely $U_{c}=[3.7,6.0,8.8]$ m/s, $F_{d,max}=0.29$ N and $\Delta
F_{d}=0.10$ N. Nevertheless, we successfully achieved a gradual variation in
$U_{c}$ while keeping the other two parameters unchanged. Note that the
optimized designs have the same rest angle and fold stiffness, but varying
radius, which is consistent with the previous results shown in Fig.2a. We
indeed demonstrated that $F_{d,max}$ and $\Delta F_{d}$ are size-independent
and, therefore, determined by the specific combination of
$(\kappa,\theta_{0})$. In Fig.4c and d, we systematically vary the jump
$\Delta F_{d}$ and peak $F_{d,max}$ respectively (target values are reported
in the Methods). Similarly, although the physical samples have drag quantities
below the target values, they still showcased the intended selective
variations. Our results confirm the validity of our inverse design approach,
but there are limitations to how closely we can match desired drag behaviors.
These limitations arise primarily from the theoretical model employed in the
optimization procedure. While it accurately captures the impact of origami
parameters on drag, it only provides semi-quantitative agreement with the
experimental data. Note that deviations between the targeted drag features and
those of the resulting optimized design in Fig.4 align with differences
reported in prior experimental-theoretical comparisons of Fig.3 (see
Supplementary Fig.7-8). Refining the model, notably using a more realistic
representation for the fluid dynamics, is thus expected to produce better
quantitative results.

##  2 Discussion

In this study, we leverage the waterbomb base’s bistability to sharply reduce
drag force with flow velocity. A uniform airflow unfolds the origami unit,
which, upon reaching a flat state, undergoes a snap-through to a more
streamlined shape, causing a sudden collapse of drag. While the use of elastic
deformation to mitigate drag increase with flow speed is not novel [1, 2, 3,
4, 5, 6], here, a decrease is achieved through significant and abrupt shape
changes over a small increment of fluid input load. Importantly, the snap-
through is tunable based on the cell radius, fold stiffness, and rest angle,
providing control over drag characteristics.
We show that the cell behavior results from the quasi-static mechanical
equilibrium between actuating fluid loading and restoring elastic forces,
captured by a Cauchy number. The snap-through occurs at a fixed critical
Cauchy number for a given initial degree of opening. The cell rest angle
serves as an additional control, influencing both snap-through onset and drag
changes. A compact cell reaches the unstable flat state at a higher $C_{y}$
and transitions to a correspondingly compact state, producing a larger drop in
drag. These experimental features are captured by a theoretical fluid-elastic
model, portraying the cell as rigid facets with folds acting as elastic
hinges, and utilizing an empirical formulation for fluid pressure forces,
consistent with prior studies. The model effectively accounts for the
influence of cell structural parameters on drag, demonstrating reasonable
predictiveness despite a tendency to overestimate drag. However, such
discrepancies are expected, considering the absence of adjustable parameters
and the simplified representation of fluid dynamics, which notably does not
consider the object’s retroaction on the flow.
The modeling framework provides the system equations to guide inverse design,
intending to identify the combination of structural parameters leading to a
targeted drag collapse. Our implementation showcases an advanced level of
control, allowing selective adjustment of key drag features: namely, the peak
drag before snapping, the drag drop, and the critical flow speed at which it
occurs. Currently, quantitative predictability is contingent on the model’s
accuracy. Improved closeness to the target is anticipated with more refined
fluid-elastic models, positioning this work as promising proof of concept. It
illustrates the potential of origami as a platform for programming drag-vs-
speed behavior.
This programmability holds value for self-protection strategies, mitigating
excessive aerodynamic loads beyond a predetermined threshold. Force
transitions above a tunable critical flow velocity could also serve as a
velocity threshold detector. The switch-like response of the waterbomb base
also holds the potential for fluidic control in a channel, for example, acting
as a relief valve. It would obstruct conduct as it deploys and discharges the
fluid beyond a predetermined pressure level, limiting pressure built up. In
all these applications, morphing is driven aeroelastically, enabling the
component to autonomously respond to the local environment, and eliminating
the need for additional sensing and control systems. Such passive mechanisms
allow for a more streamlined structural design of the aerodynamic control
surface.

Figure 5: Leveraging multistability for a double drag collapse. A multistable
origami structure is formed through the assembly of two waterbomb units, with
different fold stiffness $\kappa_{1}=19\times 10^{-3}$ N and
$\kappa_{2}=53\times 10^{-3}$ N and rest angle $\theta_{0,1}=56^{\circ}$ and
$\theta_{0,2}=71^{\circ}$, and the same radius $R=5$ cm (inset). The
successive snap-throughs lead to sequential drag collapses when exposed to a
uniform flow with increasing speed.

An advantage of origami is the ability to easily combine multiple bistable
elements to produce multistable structures [36, 37, 38, 39]. An example is
shown in Fig.5, featuring a chain made of two waterbomb cells with the same
size and different rest angle and fold stiffness (see also Supplementary
Fig.9). It produces two successive sharp transitions in the drag force. With
sufficient spatial separation, each unit is expected to have its independent
flow velocity threshold initiating snap-through. Other configurations might,
however, lead to interactions between cells through external flow, resulting
in a more complex collective response, as observed in fluidic cellular origami
(albeit for internal pressurization) [37]. Alternatively, if cells share
facets, they can also communicate the mechanical state of one unit to its
neighbor, forming logic elements [40]. This has the potential to generate more
complex snapping sequences and subsequent drag evolution with flow speed for
advanced adaptive functionalities.

## Methods

Origami fabrication. Specimens are fabricated using a layer superposition
method. The rigid facets are laser cut from a 350 $\mu$m a thick sheet of
Mylar, which is then attached to a thinner sheet using double-sided tape with
a 2.5 mm spacing that forms the creases. The folded unit is then maintained in
a desired configuration on a mount, placed in a heating chamber for an hour at
$75^{\circ}C$, and then left to cool down at room temperature. It thus
effectively prescribes the first stable state $\theta_{0}$. This annealing
process relaxes any residual stresses that might have been created during
folding [27] so that the first stable state is close to a zero-energy state.

Fold stiffness. To measure fold stiffness, a prototype of a single crease was
designed using the manufacturing process outlined above and in [30]. While
keeping one of the facets vertically fixed, we apply a variable increasing
point load to the geometrical center of the free facet and extract the opening
angle of the fold $\Psi$ (see Supplementary Fig.1a). The opening angle in the
absence of load is set to $\Psi_{0}\approx 113^{\circ}\pm 8^{\circ}$ for all
the tests, but note that the stiffness was found independent of $\Psi_{0}$ in
the previous study of [30]. The fold acts as a torsional spring, exhibiting a
linear relationship between the torque applied $T$ (accounting for the weight
of the hanging facet as well) and the angular deviation $\Psi-\Psi_{0}$ (see
Supplementary Fig.1b). The stiffness is extracted from the slope as $T=\kappa
L(\Psi-\Psi_{0})$, with $L=4$cm the crease length. Table 1 displays $\kappa$
values for the four thin sheets. Repeatability and robustness were assessed by
measuring $\kappa$ for eight folds from a 50 $\mu$m sheet, resulting in an
$8\%$ variability based on standard deviation relative to the mean value.

Sheet Thickness ($\mu m$) | $48$ | $50$ | $75$ | $100$
---|---|---|---|---
Fold stiffness $\kappa$(N) $\times 10^{-3}$ | $8$ | $19$ | $53$ | $92$
Table 1: Stiffness measured for folds made from sheets with different
thickness.

To expand the range of achievable stiffness values, we employ a strategy
inspired by prior work [34] and introduce slits in the folds. Slits with $5$
mm-length are evenly distributed along the crease length, and span its width
(see Supplementary Fig.2a). Adjusting their number modulates the effective
length of the fold, thereby changing its stiffness and facilitating the
attainment of diverse $\kappa$ values. As shown in Supplementary Fig.2b, the
stiffness of a slitted fold scales linearly with its effective length, with a
$7\%$ deviation from the expected linear relationship. To ensure that the
presence of slits has minimal impact on the interaction of the origami unit
with the flow, we conducted experiments with two specimens —one with slitted
folds and another with unslitted folds— both having the same stiffness. These
two specimens exhibited identical drag curves, indicating that the slits did
not significantly affect the performance in the flow (see Supplementary
Fig.2).

Fluid-elastic model. The equilibrium angle $\theta$ of the unit in the flow is
given by the zeros of $d\overline{(E_{el}-W)/d\theta}$, with
$\overline{E_{el}-W}$ expressed in Eq.3. As further detailed in the
Supplementary information, it yields the equation:

$\displaystyle(\psi_{m}-\psi_{m}^{0})\frac{d\psi_{m}}{d\theta}+(\psi_{v}-\psi_{v}^{0})\frac{d\psi_{v}}{d\theta}$
$\displaystyle+\frac{2}{3}C_{y}\sin^{2}\theta\sin^{2}\phi\left[A\frac{d\phi}{d\theta}\sin\theta+B\sin\phi\right]=0$
(4)

Where $C_{y}=\rho U^{2}R^{2}/\kappa$ is the same Cauchy number as defined for
experiments, $\phi$ is the angle between valley folds and the central axis,
and $A$ and $B$ are constants. Assuming rigid folding and symmetry of motion,
the expression for all angles in Eq.4 can be derived as a function on $\theta$
(see [41, 30] and Supplementary Information):

$\phi(\theta)=\begin{cases}\arccos\left(\dfrac{\cos\alpha\cos\theta}{1+\sin\alpha\sin\theta}\right)\text{
for }0\leq\theta\leq\pi/2\\\
\arccos\left(\dfrac{\cos\alpha\cos\theta}{1-\sin\alpha\sin\theta}\right)\text{
for }\pi/2<\theta\leq\pi-\alpha\\\ \end{cases}$ (5)

with $\alpha=2\pi/8$ the sector angle of a facet. Angles of the mountain and
valley folds are given by $\cos{\psi_{v}}=\cos{2\theta}$ and
$\cos{\psi_{m}}=\cos{2\phi}$. Note that the expressions as a function of
$\theta$ differ on each side of the flat configuration $\theta=\pi/2$,
describing, respectively, the opening and closing kinematics.

Eq.4 is numerically solved using the nonlinear system solver fsolve of Matlab.
The Cauchy number is varied linearly from 0 to 50 in 500 steps, and the solver
utilizes the value of $\theta$ obtained at the previous iteration as a
starting point (initializing at $\theta_{0}$ for $C_{y}=0$). As $\theta$
approaches the flat state beyond $\pi/2-\epsilon$, with $\epsilon=0.01$
radians, Eq.4 is then solved using the new set of angle expressions describing
the closing kinematics, with the first starting point at $\pi/2+\epsilon$.

From the equilibrium angle $\theta$, we obtain the dimensionless projected
surface $S/\pi R^{2}=\sin\theta\sin\phi$, and the dimensionless drag force
$\overline{F_{d}}=\pi C_{y}\sin^{3}\theta\sin^{3}\phi$ (by projecting fluid
pressure forces in the direction of the flow), for comparison with
experimental results.

Inverse design. The target values $(F_{d,max}^{t},\Delta F_{d}^{t},U_{c}^{t})$
used, along with the corresponding cell parameters, are presented in Table 2.
When varying $U_{c}$ (while keeping the other two drag features constant),
cells with the same radius but different $\kappa$ and $\theta_{0}$ are
obtained. Interestingly, selectively varying $\Delta F_{d}$ is achieved by
adjusting $\theta_{0}$ and $\kappa$ while maintaining a constant $R$.

Target values | Cell parameters
---|---
$F_{d,max}^{t}$ (N) |  $U_{c}^{t}$ (m/s) |  $\Delta F_{d}^{t}$ (N) |  $\kappa$ (N) $\times 10^{-3}$ |  $R$ (cm) |  $\theta_{0}(^{\circ})$
0.8 |  4 |  0.6 |  21 |  4.6 |  47
0.8 |  7 |  0.6 |  21 |  6.5 |  47
0.8 |  10 |  0.6 |  21 |  11.4 |  47
0.8 |  7 |  0.5 |  27 |  6.5 |  56
0.8 |  7 |  0.7 |  16 |  6.5 |  33
0.65 |  7 |  0.6 |  11 |  5.9 |  25
0.95 |  7 |  0.6 |  32 |  7.1 |  56
Table 2: Fabrication parameters produced by the optimization algorithm, for
given target drag parameters.

## References

  * [1] Vogel, S.  Drag and flexibility in sessile organisms.  _Am. Zool._ 24, 37–44 (1984).
  * [2] Alben, S., Shelley, M. & Zhang, J.  Drag reduction through self-similar bending of a flexible body.  _Nature_ 420, 479–481 (2002).
  * [3] Gosselin, F., de Langre, E. & Machado-Almeida, B. A.  Drag reduction of flexible plates by reconfiguration.  _J. Fluid Mech._ 650, 319 (2010).
  * [4] de Langre, E., Gutierrez, A. & Cossé, J.  On the scaling of drag reduction by reconfiguration in plants.  _C.R. Mécanique_ 340, 35–40 (2012).
  * [5] Schouveiler, L. & Boudaoud, A.  The rolling up of sheets in a steady flow.  _J. Fluid Mech._ 563, 71 (2006).
  * [6] Schouveiler, L. & Eloy, C.  Flow-induced draping.  _Phys. Rev. Lett._ 111, 064301 (2013).
  * [7] Wieselsberger, C. v.  Neuere feststellungen under die gesetze des flussigkeits und luftwiderstandes.  _Phys. Z._ 22, 321–328 (1921).
  * [8] Achenbach, E.  Experiments on the flow past spheres at very high reynolds numbers.  _J. Fluid Mech._ 54, 565–575 (1972).
  * [9] Lopez, D., Eloy, C., Michelin, S. & De Langre, E.  Drag reduction, from bending to pruning.  _Europhys. Lett._ 108, 48002 (2014).
  * [10] Han, J. S., Ko, J. S. & Korvink, J. G.  Structural optimization of a large-displacement electromagnetic lorentz force microactuator for optical switching applications.  _J. Micromech. Microeng._ 14, 1585 (2004).
  * [11] Tavakol, B. _et al._ Buckling of dielectric elastomeric plates for soft, electrically active microfluidic pumps.  _Soft matter_ 10, 4789–4794 (2014).
  * [12] Overvelde, J. T. B., Kloek, T., D’haen, J. J. A. & Bertoldi, K.  Amplifying the response of soft actuators by harnessing snap-through instabilities.  _Proc. Natl. Acad. Sci. U.S.A_ 112, 10863–10868 (2015).
  * [13] Shan, S. _et al._ Multistable architected materials for trapping elastic strain energy.  _Adv. Mater._ 27, 4296–4301 (2015).
  * [14] Arena, G. _et al._ Adaptive compliant structures for flow regulation.  _Proc. R. Soc. A_ 473, 20170334 (2017).
  * [15] Arena, G., Groh, R. M. J., Theunissen, R., Weaver, P. M. & Pirrera, A.  Design and testing of a passively adaptive inlet.  _Smart Mater. Struct._ 27, 085019 (2018).
  * [16] Gomez, M., Moulton, D. E. & Vella, D.  Passive control of viscous flow via elastic snap-through.  _Phys. Rev. Lett._ 119, 144502 (2017).
  * [17] Peretz, O., Mishra, A. K., Shepherd, R. F. & Gat, A. D.  Underactuated fluidic control of a continuous multistable membrane.  _Proc. Natl. Acad. Sci._ 117, 5217–5221 (2020).
  * [18] Kim, H., Zhou, Q., Kim, D. & Oh, I.-K.  Flow-induced snap-through triboelectric nanogenerator.  _Nano Energy_ 68, 104379 (2020).
  * [19] Kim, H., Lahooti, M., Kim, J. & Kim, D.  Flow-induced periodic snap-through dynamics.  _J. Fluid Mech._ 913, A52 (2021).
  * [20] Arrieta, A. F., Bilgen, O., Friswell, M. I. & Hagedorn, P.  Passive load alleviation bi-stable morphing concept.  _AIP Adv._ 2 (2012).
  * [21] Cavens, W. D. K., Chopra, A. & Arrieta, A. F.  Passive load alleviation on wind turbine blades from aeroelastically driven selectively compliant morphing.  _Wind Energy_ 24, 24–38 (2021).
  * [22] Li, S., Fang, H., Sadeghi, S., Bhovad, P. & Wang, K.-W.  Architected origami materials: How folding creates sophisticated mechanical properties.  _Adv. Mater._ 31, 1805282 (2019).
  * [23] Hanna, B. H., Lund, J. M., Lang, R. J., Magleby, S. P. & Howell, L. L.  Waterbomb base: a symmetric single-vertex bistable origami mechanism.  _Smart Mater. Struct._ 23, 094009 (2014).
  * [24] Jianguo, C., Xiaowei, D., Ya, Z., Jian, F. & Yongming, T.  Bistable behavior of the cylindrical origami structure with kresling pattern.  _J. Mech. Des._ 137, 061406 (2015).
  * [25] Yasuda, H., Chen, Z. & Yang, J.  Multitransformable leaf-out origami with bistable behavior.  _J. Mech. Robot._ 8, 031013 (2016).
  * [26] Waitukaitis, S., Menaut, R., Chen, B. G. & Van Hecke, M.  Origami multistability: From single vertices to metasheets.  _Phys. Rev. Lett._ 114, 055503 (2015).
  * [27] Lechenault, F., Thiria, B. & Adda-Bedia, M.  Mechanical response of a creased sheet.  _Phys. Rev. Lett._ 112, 244301 (2014).
  * [28] Cozmei, M. _et al._ Aerogami: Composite origami structures as active aerodynamic control.  _Composites Part B: Engineering_ 184, 107719 (2020).
  * [29] Zhang, J., Li, T., Wang, C. & Yan, X.  Aerodynamic drag characteristics of miura-ori composite structure.  _J. Aerosp. Eng._ 34, 06021004 (2021).
  * [30] Marzin, T., de Langre, E. & Ramananarivo, S.  Shape reconfiguration through origami folding sets an upper limit on drag.  _Proc. R. Soc. A_ 478, 20220592 (2022).
  * [31] Zhang, J., Changguo, W. & Zhang, L.  Deployment of smp miura-ori sheet and its application: Aerodynamic drag and rcs reduction.  _Chinese J. Aeronaut._ 35, 121–131 (2022).
  * [32] Hanna, B. H., Magleby, S. P., Lang, R. J. & Howell, L. L.  Force–deflection modeling for generalized origami waterbomb-base mechanisms.  _J. Appl. Mech._ 82 (2015).
  * [33] Overvelde, J. T. B. _et al._ A three-dimensional actuated origami-inspired transformable metamaterial with multiple degrees of freedom.  _Nat. Commun._ 7, 1–8 (2016).
  * [34] Zuliani, F., Liu, C., Paik, J. & Felton, S. M.  Minimally actuated transformation of origami machines.  _IEEE Robot. Autom. Lett._ 3, 1426–1433 (2018).
  * [35] Alben, S., Shelley, M. & Zhang, J.  How flexibility induces streamlining in a two-dimensional flow.  _Phys. Fluids_ 16, 1694–1713 (2004).
  * [36] Silverberg, J. L. _et al._ Using origami design principles to fold reprogrammable mechanical metamaterials.  _Science_ 345, 647–650 (2014).
  * [37] Li, S. & Wang, K. W.  Fluidic origami with embedded pressure dependent multi-stability: a plant inspired innovation.  _J. R. Soc. Interface_ 12, 20150639 (2015).
  * [38] Sengupta, S. & Li, S.  Harnessing the anisotropic multistability of stacked-origami mechanical metamaterials for effective modulus programming.  _J. Intell. Mater. Syst. Struct._ 29, 2933–2945 (2018).
  * [39] Jules, T., Reid, A., Daniels, K. E., Mungan, M. & Lechenault, F.  Delicate memory structure of origami switches.  _Phys. Rev. Res._ 4, 013128 (2022).
  * [40] Treml, B., Gillman, A., Buskohl, P. & Vaia, R.  Origami mechanologic.  _Proc. Natl. Acad. Sci. U.S.A._ 115, 6916–6921 (2018).
  * [41] Brunck, V., Lechenault, F., Reid, A. & Adda-Bedia, M.  Elastic theory of origami-based metamaterials.  _Phys. Rev. E_ 93, 033005 (2016).

## Data availability

The datasets generated during the current study are available from the
corresponding author upon request.

## Acknowledgments

We acknowledge support from a JCJC grant of the Agence Nationale de la
Recherche (ANR-20-CE30-0009-01) to S.R. and from the Agence Innovation Défense
to T.M.

## Author contributions

S.R., R.N. and T.M. designed the study; R.N. performed the experiments and
analyzed the results; R.N. and T.M. performed the theoretical analysis; all
authors contributed to interpreting the results and writing the paper.

## Competing interests

The authors declare no competing interests.