## Design principles

Optical imaging systems generally satisfy linear and spatial invariance, wherein the relationship between the input and output images of the system can be described through convolution operation. Specifically, for coherent light imaging systems, the input light field complex amplitude *g*in(*x*0,*y*0), output light field complex amplitude *g*out(*x*,*y*), and the system’s spatial impulse response function *h*(*x*,*y*) satisfy the equation *g*out(*x*,*y*) = *g*in(*x*0,*y*0)∗*h*(*x*,*y*), here we refer to the spatial impulse response function *h*(*x*,*y*) as a convolution kernel. Since the metasurfaces are typically 102∼103 µm in size, so we integrate it into a 4*f* system to perform the spatial spectrum modulation, of which the complex amplitude transmittance is namely the coherent transfer function (CTF), denoted by *H*(*f**x*,*f**y*), and its relationship with the convolution kernel is *H*(*f**x*,*f**y*)=Φ{*h*(*x*,*y*)}. As per the Fourier transform convolution theorem, this system can directly achieve the convolution of the input image with *h*(*x*,*y*). Consequently, diverse convolution operations can be realized by using specially designed CTF within this imaging system.

Integration and differentiation are two fundamental mathematical operations that have practical value in image processing applications for denoising and edge detection. We formulate a 3×3 binary convolution kernel depicted in Fig. 2a to perform the integral operation, whose central part, marked with a value of 0, constrains the maximum extent of the intensity smooth range. Similarly, the convolution kernel corresponding to differential operation can be expressed in another 3×3 matrix shown in Fig. 2c. Upon Fourier transformation, the normalized intensity distributions of these CTFs corresponding to two operations emerges, as illustrated in Figs. 2b and 2d. Obviously, the normalized intensity distribution of the CTF for differential operation exhibits approximate complementarity to the integral operation. To align with the demands of multichannel parallel computation, the intensity distribution of the differential operation CTF is approximated as 1−|*H*(*f**x*,*f**y*)|2 with *H*(*f**x*,*f**y*) depicting the CTF of integral operation.

In order to achieve parallel convolution operations mentioned above, we adopt the polarization encoding scheme to design a vectorial metasurface that enables wideband operation capability (*12, 14*). In the context of geometric phase metasurfaces, when the rotation angle of the polycrystalline silicon nanopillar is denoted as *θ*, for the incidence of linearly polarized light, the transmitted light induces a polarization rotation depending on the geometric phase, i.e., *ϕ* = 2*θ*. Assuming that the intensity of incident light with polarization in a certain direction is depicted as *I*0, according to the Malus' law, the intensity of the transmitted component, which coincides with the direction of incident polarization, is *I* = *I*0cos2(*ϕ*). The arrangement angle of each nanopillar unit on the metasurface can be designed according to the convolution kernel of integral operation, ensuring that the intensity transmittance of the metasurface aligns with the intensity distribution function of the integral operation CTF, i.e., *I*0cos2(*ϕ*)=|*H*(*f**x*,*f**y*)|2. For the orthogonal polarization component transmitted from the metasurface, opposite intensity distribution 1−|*H*(*f**x*,*f**y*)|2 is obtained, which approximates the CTF of differential operation. The specific design schematic of the metasurface is illustrated in Fig. 2e. From the center to the edge, the rotation angle of the nanopillar gradually varies from 0° to 45°, and the local polarization rotation of incident light field with horizontal polarization progressively changes from 0 to π/2. As depicted in Fig. 2f, this metasurface is positioned at the spatial spectrum plane of the 4*f* system. When post-selection polarization direction is consistent with the incident one, this imaging front end performs the integral operation, while implements the differential operation in the orthogonal polarization channel.

We fabricated the vectorial metasurface with transmission-type configuration by using standard electron-beam lithography and inductively coupled plasma etching. Figure 3 shows the microscope images of the metasurface, which consists of polycrystalline silicon (Poly-Si) rectangle nanopillars deposited on a glass substrate. The height and period of the nanopillar are *H* = 600 nm and *P* = 400 nm, respectively, and the refractive index is *n* = 4.018 + i0.057. The geometric size of the nanopillar is *W* = 100 nm and *L* = 188 nm, which was selected by employing the finite difference time domain (FDTD) method at 633 nm wavelength. The optical microscope image of the metasurface is shown in Figs. 3a, whose transmittance is uniform, therefore it will not affect the imaging functionality of this optical computing front end. We further characterized the transmittances of two polarization channels, whose transmissions are shown in Figs. 3c and 3d.

## Integral and differential operations

We used the experimental setup depicted in Fig. 3d to demonstrate the computational performance of the meta-device. The output light from a supercontinuum laser (SC-Pro, YSL photonics) was firstly selected by an acousto-optic tunable filter (AOTF) system, which generated monochromatic light within the range of 400 ~ 2400 nm, then passed through a polarizer (P1). The horizontally transmitted light illuminated an object (O) positioned at the front focal plane of the 4*f* system, which consists of two lenses (L1 and L2) with a focal length of 7.5 cm. The output image that contains integral and differential operation results were obtained by the CCD camera placed on the imaging plane of 4*f* system. We used the combination of a polarizer (P2) and a common COMS camera (acA2040-90 µm, Basler) to individually observe two operations. Importantly, to perform parallel computation, the results of two computing operations can be obtained simultaneously with a single exposure of a polarization camera.

As a proof of concept, the integration and differentiation were performed on the UASF1951 resolution test chart shown in Fig. 4a under 633 nm wavelength. Figures 4b and 4c show the numerically simulated results, the corresponding experimental results are shown in Figs. 4d and 4f. To facilitate a comprehensive comparison between experimental and simulated outcomes, we draw the intensity distributions along the white dashed lines of the experimental results shown in Figs. 4d and 4f, and compared it with the input and simulated results, as shown in Figs. 4e and 4g. Notably, the experimental results are consistent with the numerically simulated results, for image signal, the integral operation achieves the smooth denoising effect, while the differential operation results in edge enhancement.

From an intuitive perspective, the integral operation on an image can produce image blur and eliminate noise, and the smoothness increasing as the image size decreases. This is attributable to our designed integral operation, which assumes a value of 0 in the central region and 1 in the surrounding positions of the convolution kernel. Upon convolution operation with the incident image, the average intensity of the surrounding eight pixeled regions replaces the one of internal 0-value region, resulting in robust smoothing and noise removal effects. When the feature image size is greater than the 0-value region of the convolution kernel *h*(*x*,*y*), it does not completely remove the image blurring. Conversely, when the feature size is significantly smaller than the 0-value region, image denoising occurs. Therefore, the size of the 0-value region is critical for convolution operations on objects with a feature size. Considering the circular symmetry of the meta-imager, we take the 0-value region as a circle and define its diameter as the denoising diameter. In experiment, we set it to 88.39 µm, which is equivalent to the width of the 4th unit stripes in the second group of USAF1951 resolution test chart. Consequently, the central part of stripes with a line width > 88.39 µm retains the intensity of the original image, while stripes with a line width < 88.39 µm experience reduced or nearly complete elimination of intensity.

Performing differential operation on images will result in edge enhancement. Notably, the internal intensity of stripes in the fourth unit of the second group significantly decreases, falling below the edge intensity. For stripes with a line width > 88.39 µm, the interior is nearly completely eliminated, preserving only the edge strength. Stripes smaller than 88.39 µm experience varying degrees of intensity reduction without an evident edge detection effect. Consequently, the minimum resolution of edge detection corresponds to the denoising diameter of the integral convolution kernel mentioned earlier.

We further validated the application of these computing operations of the meta-imager in denoising and edge enhanced imaging, which were performed on specifically devised samples containing noise and onion cells at 633 nm wavelength. The denoising imaging corresponding to the integral operation was firstly tested. Figure 5a displays a flower pattern with a diameter of 1000 µm, featuring randomly distributed noise points with a diameter of 20 µm. The integral operation result is shown in Fig. 5b. Due to the denoising diameter of the meta-imager being 88.39 µm, which is much larger than 20 µm, although the noise can be completely eliminated, the edges of the flower pattern are severely blurred, making it difficult to recognize the image that should be retained. To solve this problem, a micro magnification system was incorporated into the experimental setup, enlarging the sample by 9.4 times, the denoising result is depicted in Fig. 5c. For this case, the convolution kernel diameter is smaller than the noise, the denoising effect is not significant. Further, a convolution kernel with diameter of 250 µm was used, of which the structure and corresponding CTF are shown in Figs. 5d and 5e, respectively. The integral computation result is shown in Fig. 5f, clearly, the matching convolution kernel can effectively eliminate the noise while preserving the edge contours. In practical applications, designing convolution kernel necessitates consideration not only of noise size but also other image parts and resolution of the image acquisition device. Theoretically, the greater the difference between noise and image size, the more effective the denoising achieved by the integral operation under this principle.

Finally, differentiation and integration were performed on onion cells. Figure 5g displays onion cells without computing operations. To facilitate result observable, these cells were enlarged by 9.4 times, and then differential and integral operations were carried out on these cells using a convolution kernel with diameter of 88.39 µm. Figure 5h illustrates the result of the differential operation, revealing significant elimination of the cells' interiors while preserving relatively complete edges. Figure 5i displays the outcome of the integral operation, wherein wrinkles and pits inside the cells in the original image are eliminated, and the onion cell edges are distinctly visible. The experimental results affirm that our designed meta-imager serves two operation functions: image denoising and edge detection, showing potential applications in image processing and biological sample observation.