Luyin Hu
email: luyin.hu@connect.polyu.hk

Anqing Duan
email: anqing.duan@polyu.edu.hk

Mengying Li

David Navarro-Alarcon

Andrea Cherubini
email: andrea.cherubini@lirmm.fr

Lu Li

i Paint with the Sun: A Thermal-Vision Guided Robot to Harness Solar Energy for Heliography Luyin Hu, Anqing Duan, Mengying Li, Andrea Cherubini, Lu Li, and David Navarro-Alarcon Abstract-In this paper, we present a novel robotic system for heliography, which literally refers to painting with the sun. In the context of artistic creation, heliography involves the precise manipulation of a magnifying glass to burn artwork on wooden panels by concentrating solar energy. Since this painting procedure exposes the human operator to prolonged periods of intensive sunlight, it can potentially cause sunburns to the skin or damage the eyes. To avoid these issues, in this paper we propose to automate the heliography process with a robot in lieu of a human operator. The proposed thermal servoing capabilities have the potential to robotize various solar power generation technologies such as Concentrated Solar Power (CSP) and Concentrated Photovoltaics (CPV). To perform this task, our robotic system is equipped with a magnifying glass attached to the end-effector, and is instrumented with vision and thermal sensors; The proposed sensor-based controls enable to automatically: 1) Track the orientation of the sun in real-time to maximize the concentrated solar energy; 2) Direct the solar rays towards the point of interest; 3) Control the heat power intensity at the target point to achieve the required steady state temperature. The performance of our system is evaluated by conducting autonomous heliography experiments with various patterns. Index Terms-Thermoception; Visual servoing; Sensor-based control; Robotic manipulation; Solar energy.

I. INTRODUCTION

I N the context of art production, heliography is a technique that consists in concentrating-with the aid of mirrors and lenses-solar rays on wooden surfaces to create artworks [START_REF] Kruse | Sunlight art: Painting with a magnifying glass[END_REF]. Unlike painters, who draw on a canvas using a paintbrush, heliographic artists rely on the focal point that results from the refraction of sunlight passing through a magnifying glass. This focal point possesses a high energy flux that enables to carbonize the wood and leave patterns, similar to a paintbrush. Fig. 1(a) shows an example of a traditional heliographic artwork produced by manipulating the focal point on a wood canvas. Such a procedure requires precise manipulation skills, since both the focus duration and the projection distance affect the quality of the artwork; A long focus duration will burn through the wood, while a deviated projection distance will not produce the desired width of the burning trace. To produce a complete art piece may take an artist anything between a couple of hours to several weeks, depending on the size, level of detail, and employed materials. While creating the desired pattern, the artist may undergo prolonged periods of direct exposure to sunlight and its ultraviolet radiation, which can cause various skin and eye problems, e.g., sunburns, melanoma, or even cataracts. As a feasible solution to these issues, we propose to automate heliography with a robotic system that controls the sun concentration process, effectively taking over the role of the artist, see Fig. 1(b). The sensorimotor controls that are needed to produce this type of artworks are also valuable from a renewable energy perspective, e.g., robotic structures can automate the harnessing of solar energy in a concentrated solar power system. Conceivably, robots with the capability to exploit off-the-peg energy sources from the environment can potentially save a considerable amount of energy. Our aim in this paper is precisely to develop the necessary models and methods for robots to achieve this level of thermomotor intelligence.

From an application perspective, there is a large volume of literature dedicated to robotizing tasks such as painting, drawing, and brushing, whose spatial trajectories are similar to those used in heliography. For example, to teach a robot the skill of performing aesthetic calligraphy, a controller can reproduce the strokes captured from previous human writing demonstrations [START_REF] Sun | A calligraphy robot-callibot: Design, analysis and applications[END_REF]. Similarly, artistic portraits can be realized by processing images of human faces and then planning trajectories to be followed by a robot-held drawing device [START_REF] Calinon | A humanoid robot drawing human portraits[END_REF]. Different from the aforementioned painting robots, our system does not need any ink or pigment to leave traces on a canvas, instead, it relies on the carbonization of a wooden panel. Similar working principles can be found in systems that manipulate lasers [START_REF] Muddassir | Robotics meets cosmetic dermatology: Development of a novel visionguided system for skin photo-rejuvenation[END_REF], which have been used to create various types of art works [START_REF] Tamrin | Laser discoloration in acrylic painting of visual art: Experiment and modeling[END_REF], [START_REF] Xyntarianos-Tsiropinas | Photon painter: A robotic device that transfers drawings onto a large surface area using laser beam on photomechanical reaction materials[END_REF]. Compared with these laserbased painting systems, our robot does not require any external power supply to carbonize the surface.

We highlight that one salient aspect of our proposed system is exploitation of solar energy, an important topic in various fields due to its sustainable and eco-friendly nature. One popular application is the so-called solar cooker that can process food or boil water by converting solar energy into heat [START_REF] Yettou | Solar cooker realizations in actual use: An overview[END_REF]. Our system exploits solar energy to carbonize the wooden panel, and unlike the solar cooker, it needs faster and more precise temperature controls to avoid burning through the wood. Many mobile robotic systems power their motion by adjusting the pose of onboard photovoltaic panels so as to maximize the incidence of solar rays over its surface [START_REF] Plonski | Energy-efficient path planning for solar-powered mobile robots[END_REF], [START_REF] Gorjian | Applications of solar pv systems in agricultural automation and robotics[END_REF], a condition analogous to our rays-lens alignment problem. The heliography process requires to maximize the amount of solar rays passing through the magnifying glass, which is achieved by tracking the sun's direction [START_REF] Chu | Sun-tracking imaging system for intra-hour dni forecasts[END_REF]. To this end, we design an iterative strategy that estimates the sun's direction by visually measuring the projected shadow feature of the lens; At the core of this approach is the system's interaction matrix that captures the relationship between the sun direction increments and the feature vector increments. Compared with alternative approaches for sun direction inference based on specialized sensors (digital compass [START_REF] Afarulrazi | Solar tracker robot using microcontroller[END_REF], catadioptric camera [START_REF] Kadmiri | A novel solar tracker based on omnidirectional computer vision[END_REF], light dependent resistors [START_REF] Arlikar | Three dimensional solar tracker with unique sensor arrangement[END_REF]), or based on deep neural networks (which require large datasets [START_REF] Peretroukhin | Inferring sun direction to improve visual odometry: A deep learning approach[END_REF]), our proposed sensor-based solution is simple to implement and data efficient.

For the robot to automate the task, real-time acquisition of the surrounding environment states is essential [START_REF] Ma | Application of sensor data information cognitive computing algorithm in adaptive control of wheeled robot[END_REF]. In our application, two states are considered to play an important role in the success of the task: i) location of the solar focal point, and ii) temperature of the projected focal point over the surface. The former suggests installation of an RGB camera to visually track the projected 2D point [START_REF] Cippitelli | Radar and rgb-depth sensors for fall detection: A review[END_REF] while the latter suggests installation of a thermal camera to monitor its temperature [START_REF] Shor | Compact thermal sensors for dense cpu thermal monitoring and regulation: A review[END_REF], [START_REF] Hu | On radiationbased thermal servoing: New models, controls and experiments[END_REF]. The integration of thermal and visual information to augment robot capabilities has been recently studied in many works, e.g., human falling detection [START_REF] Zhong | Multi-occupancy fall detection using non-invasive thermal vision sensor[END_REF], thermographic reconstruction [START_REF] Cao | Depth and thermal sensor fusion to enhance 3d thermographic reconstruction[END_REF], spacecraft operations [START_REF] Palmerini | Combining thermal and visual imaging in spacecraft proximity operations[END_REF] etc. In this work, we show how thermal and visual signals can be used to guide the end-effector's motion, thus, paving a new path towards the development of multi-perception modalities for robotic manipulation [START_REF] Cherubini | Sensor-based control for collaborative robots: Fundamentals, challenges, and opportunities[END_REF].

The original contributions of this paper are as follows:

• A new sensor-based sun tracking algorithm to collect a maximum amount of solar energy with the lens; • A new motion planning strategy to direct and concentrate the sun's rays over a point of interest; • A new sensor-based control law to regulate the temperate at the target point and achieve a desired steady state or an appropriate heating process.

II. PRELIMINARIES Notation Throughout this manuscript, we denote all column vectors by lower bold case letters, e.g. r ∈ R n×1 , and matrices by bold capital letters, e.g. T ∈ R m×n . We use a left superscript to indicate the coordinate system of a homogeneous position vector, e.g. a r ∈ R 4×1 is expressed in frame (xyz) a . The axes of (xyz) a are denoted by x a , y a , z a , respectively. A matrix T a b ∈ R 4×4 describes the homogeneous transformation from frame (xyz) a to frame (xyz) b . The trigonometric functions sin θ and cos θ are abbreviated as Sθ and Cθ, respectively.

A. Problem Formulation

Fig. 2 provides an overview of the proposed workflow of our method. The names of the blocks are italicized in the rest of this section. The robotic heliography process consists in concentrating solar energy at a group of target positions in the world frame P w = { w p * 1 , . . . , w p * k }. These positions are obtained by performing feature extraction and automatic scaling on an input image I t (viz. the target figure), which generates the corresponding ordered group of points [START_REF] Huo | Keypoint-based bimanual shaping of deformable linear objects under environmental constraints using hierarchical action planning[END_REF]. Humans achieve this task by manipulating the lens' configuration according to past experience and current observations. To reproduce such functionality with a robot, we analyze the lens' optical principle and developed a ray tracing optical simulation. This simulation sets the ground for an imagebased method to estimate the solar orientation and to compute the required lens' configuration that concentrates solar energy at w p * i . To create the desired patterns, we develop a path planning algorithm which takes into account the operation time and the overlap between subsequent target positions.

Alongside the coupling between the robot configuration and the position where solar energy is concentrated, another key aspect of our problem is the heat transfer process occurring at w p * i . The charred regions that appear on the wood's surface result from pyrolysis [START_REF] Rath | Heat of wood pyrolysis[END_REF], a reaction that is triggered when the temperature exceeds a certain value. In heliography, it is critical to track the pyrolysis process to avoid under/overheating the surface; Humans detect this by observing the emerging smoke. However, it is complicated for a robot to visually detect smoke due to the intense light arising from the process. Instead, in our proposed method we equip the robot with a near infrared thermal camera to directly measure the surface's temperature profile. This sensor feedback is processed to search for the maximum temperature value, which in then sent to the robot's motion controller.

B. Proposed Sensor-Based Robotic System

The main components of the proposed system are depicted in Fig. 3. Two servo motors are attached to the end-effector of a 3 degrees of freedom (DOF) robot arm Dobot Magician, to enable the control of 5-DOF poses. A spherical Fresnel lens is manipulated by the robot to concentrate solar energy, with an RGB camera and a thermal camera rigidly attached to it to monitor the process. The mechanical structure of the end-effector is designed so that the optical axis of the Fresnel lens and the optical axes of the two cameras intersect at the theoretical focal point. The lens' yaw and pitch angles are independently controlled by motor 1 and motor 2, while its 3-DOF translation is controlled by the robotic arm. The joint configuration of the system is denoted by x ∈ R 5 . Based on the robot's kinematics, we can set the distance between the lens' center and a target point over a plane along with the lens' axis 1 .

We denote the static robot base frame by (xyz) r and the world frame (which is determined by an ArUco marker attached to the ground plane) by (xyz) w . We denote the origin of the lens frame (xyz) l by O l and set it at the center of the lens upper surface. The axes of (xyz) l are aligned with the axes of (xyz) r , which defines the origin of the lens yaw angle (the angle between the axis y l and the plane x r y r ) and the lens pitch angle (rotation around the axis x r ). The RGB camera frame (xyz) c1 and the thermal camera frame (xyz) c2 are set according to the conventions from [START_REF] Andrew | Multiple view geometry in computer vision[END_REF]. Based on the 1 Provided that such desired pose is within the robot reachable workspace. system's kinematics, the robot-lens transformation matrix T r l can be derived. Since the cameras and the Fresnel lens are rigidly attached to the end-effector, T l c1 , T l c2 , T c1 c2 are known constant matrices. At the beginning of each solar painting task, the world-robot transformation matrix T w r is obtained by observing the ArUco marker attached to the ground plane.

III. SOLAR ORIENTATION ESTIMATION

A. Geometric Optics

In this section, we present the development of our ray tracing model, which is a critical part of the motion controller. We include the necessary definitions and derivations to simulate ray tracing in solar painting tasks. Interested readers may refer to [START_REF] Lin | Advanced geometrical optics[END_REF] for more information.

The manipulated Fresnel lens relies on the same working principle as traditional spherical convex lenses. Yet, it is more efficient in terms of cost and energy transmission, ascribed to its specialized design [START_REF] Xu | Simple fiber-optic refractive index sensor based on fresnel reflection and optical switch[END_REF]. As depicted in Fig. 4(a), the curvature design of a convex lens is based on a sphere with center O o and radius r s as reference. The region of interest (colored in turquoise) is confined by the lens azimuthal angle θ l ∈ [0, 2π) and by the lens polar angle ϕ l ∈ [0, ϕ lb ], where ϕ lb is a constant lens bounding polar angle. Fig. 4(b) shows that the incident and refracted rays are denoted by o s 0 , o s 1 , o s 2 under the optics frame (xyz) o . The rays only refract at the boundary surface, hence, the inner part (colored in light blue in Fig. 4(c)) barely contributes to concentration, thus, it can be removed to alleviate energy loss. The final design of the Fresnel lens exploiting this feature is shown in Fig. 4(d), where the critical curved surface (colored in green) is divided into pieces and placed on the flat boundary surface.

Following the conventions in [START_REF] Lin | Advanced geometrical optics[END_REF], the origin of the optics frame (xyz) o is set at the center of the reference sphere. We set the axes of (xyz) o to be parallel to the axes of (xyz) l . Since the thickness of the Fresnel lens is small (2 mm in this study), to simplify the derivation we fairly assume that the sun ray enters and leaves the lens at the same point (See Fig. 4(d). Consequently, the position of an incident point o p l1 and its corresponding leaving point o p l2 are expressed (using homogeneous 4-vectors [START_REF] Andrew | Multiple view geometry in computer vision[END_REF]) as follows:

o p l1 = o p l2 = p l1x p l1y p l1z 1 = r s Sϕ l Cθ l r s Sϕ l Sθ l r s Cϕ lb 1 . (1)
The directional incident ray is denoted as:

o s 0 = Sϕ 0 Cθ 0 Sϕ 0 Sθ 0 Cϕ 0 0 , (2)
where ϕ 0 ∈ [ π 2 , π) and θ 0 ∈ [0, 2π] are the ray polar angle and the ray azimuthal angle (see Fig. 4

o n l1 = Sϕ l Cθ l Sϕ l Sθ l Cϕ l 0 , o n l2 = 0 0 1 0 . (3)
The incident angle between o s 0 and the convex boundary normal o n 1 is calculated as

γ 1 = arccos(-o s 0 • o n l1 ).
To obtain the refracted ray o s 1 , we first define:

ω 1 = N 1 Cγ 1 -1 -N 2 1 + (N 1 Cγ 1 ) 2 (4)
for N 1 = ξ 0 /ξ 1 with ξ 0 and ξ 1 as the refractive indices (relative to a specific light wavelength) of the air and of the lens material, respectively. Then, o s 1 can be calculated as [START_REF] Lin | Advanced geometrical optics[END_REF]:

o s 1 = s 1x s 1y s 1z 0 = ω 1 o n l1 + N 1 o s 0 (5)
while the refracted ray o s 2 can be similarly calculated as in the procedure above. After the ray leaves the lens's flat boundary, it travels rectilinearly until hitting the ground plane, whose normal direction is denoted by o n g . We denote the intersection point of o s 2 with this plane by:

o p g = o p l1 + K o s 2 (6)
where K represents the scaling factor of o s 2 . For a point o p go = [0, 0, d, 1] in the ground plane, we can calculate K by solving the equation

( o p g -o p go ) • o n g = 0.
This enables to obtain the analytical expression for o p g , which is a function of the following variables: Remember to change r l to r l i

1) The direction of the incident ray o s 0 , which depends on ϕ 0 and θ 0 according to [START_REF] Sunscribes | [END_REF]. In this study, ϕ 0 is related to the sun elevation angle and θ 0 is related to the sun azimuthal angle, which vary with time. 22) The point where the ray hits the lens, which is parameterized by ϕ l and θ l . In the next section, we introduce a group of incident positions that simplify its computation.

3) The relative configuration between the lens plane and the ground plane, which is described by T c1 l and can be controlled by the robot's motion.

B. Ray Tracing Simulation Method

To analyze and exploit the coupling between the controllable lens configuration and the solar orientation, we developed a visualization tool based on the aforementioned derivation with the Python Plotly library [START_REF] Inc | Collaborative data science[END_REF]. An illustration of our ray tracing method is shown in Fig. 5, where the solar incident angles are set to ϕ 0 = 5π/6, θ 0 = 0, the ground plane normal is set to o n g = [0, 0, 1, 0] , and the lens surface is parallel to the ground plane. To generate a set of uniformly distributed incident points, we pick a series of lengths r li that are evenly distributed from 0 to r s Cϕ lb (the lens radius), where different values of r li are represented by different colors. For each circle of radius r li , we select sampling points with an equal adjacent θ l to be visualized. The values of θ l of the selected points are shown with varying color saturations. As shown in Fig. 5(a), the parallel sun rays first converge then diverge after passing through the lens. Fig. 5(b) depicts the incident points on the circular lens surface, while Fig. 5(c) depicts the projected points over the ground plane.

Our ray tracing method shows that for a group of incident points o p l1 with the same r li , their corresponding projected points o p g form a (non-circular) closed contour over the ground plane. Note that for the group of incident points with the largest r l (depicted in red in Fig. 5(b)), the corresponding projected points also form the outermost closed contour in Fig. 5(c). These correspondences allow us to only analyze a representative group of incident points with the same r li . Based on this observation, we utilize the geometric characteristics (centroid and orientation) of the outermost boundary of the projected light spot for solar orientation estimation, which can be computed from the visual feedback.

C. Image-Based Solar Orientation Estimation

To automatically align the lens with the sun rays direction, we developed an efficient method that estimates the sun's orientation based on feedback images. To this end, it is important to note that the ground projection point o p g is coupled with the unknown sun orientation and the known lens-ground relative configuration. Our strategy is to estimate the solar orientation based on images of this point, and on the kinematic relationship between o p g and ϕ 0 , θ 0 , which is obtained from the ray tracing model.

Consider a 2D image captured by the camera mounted on the end-effector. We segment the image pixels that represent the projected sun spot by first transforming the image into gray scale and then applying a binary threshold to find the spot's region [START_REF] Huo | A sensor-based robotic line scan system with adaptive roi for inspection of defects over convex free-form specular surfaces[END_REF]. To relate the real-world sensor data with the numerical model, we must compute o p * g , i.e., the 3D coordinates in the optics frame (xyz) o of these image pixels.

Let i p ∈ R 3 be the homogeneous coordinates of a relevant 2D pixel point, and w p = [x w , y w , z w , 1] its corresponding unknown 3D point expressed homogeneous coordinates. According to the pinhole camera model [START_REF] Zhang | A flexible new technique for camera calibration[END_REF], these two points satisfy the following relation:

i ps = A r 1 r 2 r 3 t w p, (7)
with s = 0 as a scaling factor, A ∈ R 3×3 as the intrinsic parameter's matrix, and T w c1 = [r 1 , r 2 , r 3 , t] ∈ R 3×4 as the extrinsic parameter's matrix. Since the projected points are in the same plane (viz. the ground plane), i p and w p satisfy the following homography transformation:

x w y w 1 = A r 1 r 2 t -1 i ps. (8)
This way, we can compute o p img , the corresponding point of an image pixel i p under the optics frame via:

o p img = T o l T l c1 T c1 w x w y w 0 1 . ( 9
)
Once the image pixel has been transformed to the frame (xyz) o , we can formulate a simulation-based servoing problem to estimate the unknown solar orientation. For that, we select the outermost boundary as a representative image contour. We denote the center of this ellipse on the ground plane by o c img = [x c , y c ] , from which the feature vector y ∈ R 2 of polar coordinates is computed as:

y = Π Θ = x 2 c + y 2 c 1 2 arctan(y c /x c ) (10)
The construction of y is based on empirical observations of the ray tracing model. In its general form, the elements of y is coupled with solar polar angle ϕ 0 and solar azimuthal angle θ 0 , such that y = y(Π(ϕ 0 , θ 0 ), Θ(ϕ 0 , θ 0 )). When the lens plane is parallel to the ground plane, the variation of ϕ 0 affects the distance Π, and the variation of θ 0 affects the angle Θ, see Fig. 6, such that y = y(Π(ϕ 0 ), Θ(θ 0 )).

The estimation of the solar orientation computed at the time step t is denoted by the vector ϑ t = [ϕ t 0 , θ t 0 ] , and its corresponding feature vector by y t = [Π t , Θ t ] ; These two structures satisfy y t = y t (ϑ t ). To differentiate the feature vector computed from image feedback with the feature vector computed from simulation, we denote the former by y img and the latter by y t . Our method estimates the solar orientation by iteratively minimizing the difference between y t and the feature y img computed from image measurements (using ( 10)) as follows:

ϑ t+1 = ϑ t -λA -1 t (y t -y img ), (11)
where λ > 0 is a gain, and A t ∈ R 2×2 is the interaction matrix of the differential relation d dt y t = A t d dt ϑ t , and is numerically computed with the ray tracing model as follows:

A t =    Π t (ϕ t 0 + δ, θ t 0 ) -Π t δ Π t (ϕ t 0 , θ t 0 + δ) -Π t δ Θ t (ϕ t 0 + δ, θ t 0 ) -Θ t δ Θ t (ϕ t 0 , θ t 0 + δ) -Θ t δ    (12)
for δ > 0 as a small constant. The interaction matrix can be further simplified, e.g., if the RGB image is acquired when the lens surface is parallel to the ground plane, the variation of θ t 0 has no affect on Π t , and the variation of θ t 0 has no affect on Θ t , thus, A t is diagonal. Also, note that we should set an initial guess for ϕ t 0 = 0, to avoid singularities of A t at start. The estimation algorithm terminates when y ty img is smaller than predefined value. The final estimation at the time instance t f of the solar angles is denoted by ϑ = [ ϕ 0 , θ 0 ] = ϑ t f .

The selection of the geometric features for solar orientation estimation is not unique. There are other possible alternatives, e.g., the orientation or the shape of the ellipse. Since the solar orientation is parameterized by two variables, utilizing two features for estimation is sufficient. One requirement of the selected features is to be independent, to avoid A t from becoming singular. The selected features should also be robust enough, so that the estimation is reliable even in the presence of noise and inaccuracies, as we will discuss in Section V-B. We select features Π and Θ for our algorithm, based on the aforementioned criteria. .

IV. MOTION CONTROL ALGORITHMS

A. Lens Configuration Control

In this section, we derive a method to control the focus of solar rays over the target point w p * i . To this end, we model the solar energy power incident on the lens surface as:

w l = q s A l cosβ dA l , (13)
where q s is the solar irradiance (power per unit area), A l the area of the lens surface, and β the angle between the sun ray and the lens surface normal. Note that as the lens is flat, β remains constant for every dA l . To maximize the solar energy power incident to the lens, the lens surface normal should be aligned with the estimated orientation of the incident sunlight.

The estimated solar orientation in the robot frame (xyz) r is calculated as r s 0 = T r l T l o o s 0 , where o s 0 is the estimated solar orientation computed by substituting ϑ into (2). The target point w p * i is the location on the ground where the concentrated light spot needs to be placed, which is depicted in Fig. 3; It can be transformed into the robot frame (xyz) r as r p * i = T r l l r c1 T c1 w w p * i . The target lens position in the robot frame (xyz) r is then computed as:

r p * l = r p * i -d lp r s 0 , (14)
where d lp is the distance between the lens center and the target point. The following equation describes how to align the orientation of the lens towards the estimated solar orientation r s 0 , which is expressed in zyx Tait-Bryan angles α * l , β * l , 0:

-Cα * l Sβ * l -Sα * l Sβ * l Cβ * l 0 = -r s 0 (15)
The target lens pose is constructed as

Ω * = [ r p * l , α * l , β * l , 0]
. By using the robot's inverse kinematics we can then transform the target pose Ω * into its corresponding target joint configuration x * , which is then sent to the robot's control system to command the motion. Fig. 7 depicts the effect of varying the lens configuration on the concentration of sun rays. This visualization is divided in various sequential processes, i.e., the alignment of the polar angle, the alignment of the azimuth angle, and the translation of the lens to r p * l . This figure shows that by first orienting the lens parallel to the ground and then translating it towards the incident point, the collected solar power and the lens' optical concentration can both be optimized. This strategy is adopted by our proposed motion controller.

Nevertheless, we argue that it is not strictly required to include the reorientation process for all robot solar concentration tasks. As shown in Fig. 7 (b), the lens is capable (to some extent) of concentrating solar energy even if the lens is not perfectly parallel to the ground. When the robot's DOF are limited, it is still possible to achieve suboptimal solar concentration. Moreover, from Fig. 6(a), we can see that the distance between the lens and the optimal focal point decreases while the angular offset between lens surface normal and solar elevation increases. This unique optical characteristic could be exploited by robots with limited manipulability.

B. Parameter Selection for Motion Control

In this section, we analyze the heat transfer process involved before, during, and after wood pyrolysis. Then, we design a series of experiments, to identify the motion control parameters, specifically, the desired lens configuration and the termination criteria for the solar pyrolysis process.

1) Pre-Pyrolysis Stage: Wood pyrolysis is a chemical reaction that is triggered when the wood temperature T w reaches the pyrolysis point T p . Depending on the type of wood, T p varies from 200 °C to 400 °C. During the heating process, if T w < T p , the external heat inflow mainly induces physical changes, including the vaporization of the volatiles, and the increment of the wood's temperature T w . According to the heat transfer principles, if the external heat flux q e (power per unit area) is not strong enough, T w will only reach a steady state temperature T w = T ss < T p , but no charred regions will be produced. To actually trigger the desired reaction at the target point, q c should be maintained above a target threshold q * c . In our study, the relationship between q c and the robot joint configuration x is modeled as:

q c = ακq s A l /A c (x), (16)
where α ∈ (0, 1) denotes the absorptance of the ground plane, κ ∈ (0, 1) a factor quantifying the solar energy loss due to lens absorption, A l the lens surface area, A c the concentrated solar spot area, and A c (x) the function mapping the robot joint configuration x to A c . Since the lens configuration is always aligned with the solar orientation, we can vary q c by changing d lp (i.e., the distance between lens center and target point), so that A c (x) ≡ A c (d lp ). The mapping function A c (d lp ) is calculated with our proposed ray tracing simulation model. However, note that it is difficult to accurately obtain the true parameters α, κ, q s , q * c . To deal with this issue, instead of relying on adaptive temperature control [START_REF] Hu | On radiationbased thermal servoing: New models, controls and experiments[END_REF], we implement a simple yet effective solution where we set d lp = d lp , for d lp as an empirical value, based on collected experimental data. To this end, we conduct a series of experiments with varying d lp and a fixed heating duration of around 20 seconds. With our tests, we found that the target point is carbonized when d lp is set to 8.2 cm and 8.6 cm, which is close to the focal length indicated by the manufacturer (8.0 cm). To terminate the heating process, we set d lp = d lp + δ d lp where the offset is set to δ d lp = 3 cm. The temperature profile at the target point is monitored with the thermal camera, and its maximum value at time instant t is denoted by T max t . Fig. 8(a) depicts the temporal evolution of this variable. Note that in these tests, T max t is capped at T max c = 180.9 °C, which is the hardware upper limit of our thermal camera. Since T max c is close to the theoretical range of wood pyrolysis temperature, we can fairly assume that T max c = T p to compensate for the hardware's limitations.

2) Pyrolysis Stage: Once pyrolysis has been triggered, the wood components (e.g., cellulose) start to react with oxygen to produce water, carbon dioxide, and the charred spots. Afterwards, if q c is continuously applied, the pyrolysis region starts to expand over the surface due to the exothermic chain reaction and propagation of vaporization. In Fig. 8(b), we depict successive thermal images collected during the tests with d lp = 8.6 cm. We find that the temperature of the target point's center region increases faster than its surrounding, and it is the first to reach T p . We denote by N p the number of pixels having reached T p in a thermal image, and depict the temporal evolution of this area for d lp = 8.6 cm in Fig. 8(c). We choose N p as a proxy of the actual carbonized area of the target point.

As shown in Fig. 8(a), it is more time-efficient to create a relatively small carbonized area by setting d lp = 8.2 cm, than by setting it to d lp = 8.6 cm. To create a series of carbonized points of uniform size, in the motion control algorithm we set N p ≥ N p as the termination criteria for the solar concentration method. Our tests reveal that setting N p = 30 will produce a carbonized point with a diameter of around 0.2 cm. 3) Post-Pyrolysis Stage: In this study, the surface area of the Fresnel lens is relatively small, therefore, the concentrated solar energy inflow is not strong enough to trigger a spontaneous exothermic process (i.e., the continuous burning of the wood). Thus, when the process of the solar concentration is terminated by the robot motion, the temperature of the target point should drop immediately. However, we find that for tests where d lp = 8.2 cm and d lp = 8.6 cm, T max t remains at a constant high value for around 3 seconds after the solar concentration is terminated. This is mainly due to the actual pyrolysis temperature being higher than the camera's upper limit T max c , as it takes time for T w to drop to T max c ; Afterwards, the actual variation of T w can be captured by the thermal camera. This phenomenon is worth mentioning, since we are using N p as the termination criterion of the heating controller. In such case, if we start to heat the next target point immediately after the condition N p ≥ N p is fulfilled, the high temperature region of the previous target point will also be counted in the current calculation of N p . To solve this problem, we programmed the robot to not proceed to heat the next target point until N p has dropped to null.

C. Heliography Path Generator

The aim of the developed robotic system is to create uniform carbonized black points at a group of target positions P w on the ground plane. To this end, the first step is to transform the target image into a binary image I b . Since the resolution of I b depends on I t , if each pixel with a positive value in I b is mapped to P w , the operation time of a solar painting task will be very long, and some target positions might be out of the robot's working range. Moreover, since the solar energy is concentrated at the area of interest A c , we require an algorithm for avoiding overlapping elements in P w . To solve the aforementioned problems, feature extraction and scaling operations are applied to I B to generate a low-resolution feature image I f . Feature extraction techniques are applied depending on the characteristics of I t . For monochrome target images including text, logos, and symbols as in Fig. 9 (a 1 , b 1 , and c 1 ), we use Canny edge detector from the OpenCV library [START_REF] Bradski | The OpenCV Library[END_REF] to extract the contours in the image as a representative feature. Note that in Fig. 9 c 1 , a part of the trace is so narrow that it may cause path overlapping. Therefore, we use the center-line of the closed contour as desired feature. For target images that contain abundant color information like the humanoid portrait shown in Fig. 9 (d 1 ), the pattern of the images cannot be well represented by edges; In this case, we use the adaptive threshold function and a series of denoising filters implemented in OpenCV, to generate I b . Now, we introduce the method for scaling and mapping the image pixels to the pattern's points. Consider the binary image I b of dimensions x b ×y b (pixels). We define the working canvas on the ground plane as a square region of dimensions l w × l w cm. The scaling factor is computed as:

ω s = l w /(max(x b , y b ) • r p ) (17)
where r p is the radius of the projected light spot. As mentioned in Sec. IV-B, the parameter of the motion controller is set to N p = 30, which will produce a carbonized point with a diameter of around 0.2 cm. Thus, we set r p = 0.2 cm to compute the scaling factor ω s . For each image pixel (u, v) from I b , its corresponding target position on the ground plane is computed as w p * i = [ur p + d x , vr p + d y , 0, 1] , where (d x , d y ) defines the displacement from the origin (the center of the ArUco marker). These computed world coordinates are then grouped according to u, ordered according to v and stored in an ordered list P w , which is then sent to the motion controller described in Algorithm 1.

Recall that in Sec III-C, the solar orientation with respect to the robot frame is estimated and denoted by ϑ. As the actual solar configuration varies with time during a solar painting process, ϑ should be updated accordingly. In our method, instead of solving the simulation-based servoing problem again (as it requires additional computational time), we update ϑ by adding a time-varying correction factor. To this end, let us denote the the conventional solar orientation (i.e., the azimuthal/elevation angles with respect to north/horizon) at the time instant t by ϑ t c . Given the latitude and longitude information, ϑ t c can be obtained from the Python library Pysolar [START_REF]Staring directly at the sun since 2007[END_REF]. At the system calibration stage, the initial conventional solar orientation is recorded as ϑ t0 c . Accordingly, the solar orientation with respect to the robot frame at time instant t, the orientation ϑ can be updated as:

ϑ t = ϑ + (ϑ t c -ϑ t0 c ), (18)
where the term ϑ t c -ϑ t0 c effectively acts as the time-varying correction factor.

V. RESULTS

A. Focal Length Calibration

In this section, we compare the results of A c computed from simulation with those from the sensor feedback. To conduct the validation, we first align the orientation of the lens with the estimated sun orientation, then, control the robot to focus the sun on a target point with varying lens-plane distance d lp from 6.4 cm to 9.8 cm, with a uniform step of 0.1 cm. As shown in Fig. 10(a), the contour of the projected light spot is extracted and depicted in red. The image pixels of the solar spot are then transformed into 3-D coordinates in (xyz) r coordinates, and the actual area of the solar spot is computed accordingly. Note [START_REF] Shor | Compact thermal sensors for dense cpu thermal monitoring and regulation: A review[END_REF]; end that this step is necessary, since the position of the camera is also moving with respect to the ground plane, making the direct correlation between the amount of the pixels of interest in the image and the actual solar spot size ambiguous.

i ← i + 1; Update ϑ t ← (
The relationship between the actual A c computed from experiments and the lens-plane distance d lp is depicted in Fig. 10(c) with red dots. The theoretical values of A c under the same solar orientation and lens configuration are obtained through the simulation and depicted in Fig. 10(c), with the blue dashed line. We find that the coupling between A c and d lp is similar in the experiment and the simulation, and the value of d lp that corresponds to the minimum A c in the experiment is 8.2 cm, which is close to the focal length f l = 8.0 cm provided by the manufacturer. The match between theory and practice demonstrates the effectiveness of the proposed sensorbased calibration method. Nevertheless, obtaining A c from real-time sensor feedback is a challenging task from an image processing perspective. As shown in Fig. 10, there may exist regions of the image (e.g., the white region of the ArUco marker) with similar brightness as the concentrated light spot.

B. Solar Orientation Estimation Experiments

We conducted four experiments denoted by E i=1,2,3,4 . In each experiment, the orientation of the robot is arbitrarily set; The lens plane is set to be parallel to the ground plane, the lens yaw angle is set to vary within a range β * l ∈ [-B i , B i ] with a step of 0.5°, where B 1 = B3 = 25°, B 2 = 35°, B4 = 15°, while the lens center O l is maintained at a constant position. In these estimation experiments, we collected around 400 RGB images of the light spot. The k th captured image in the E i experiment is denoted by I i k , and the corresponding robot configuration x i k is also recorded. Some representative images I i=1,2,3,4 k are depicted in Fig. 11 (a 1 )-(a 4 ), respectively, and the corresponding robot configurations are visualized in Fig. 11 (c 1 )-(c 4 ), respectively. Fig 11 (b) provides the annotations of the robot configuration in these tests. The term "original end-effector" refers to the tool center point of the robot arm where the motorized lens system is added, which is denoted as O oe in Fig. 3.

In Fig. 11 (b) and (c 1 )-(c 4 ), the red ellipse representing the light spot is generated by mapping the contour of interest in I k i to the robot frame (xyz) r . The (position-based) geometric feature vector y of the red ellipse, which is defined in Eq. [START_REF] Gorjian | Applications of solar pv systems in agricultural automation and robotics[END_REF], is used to estimate the solar azimuthal angle θ 0 and the solar elevation angle ϕ 0 3 . To quantitatively evaluate the performance of the solar orientation estimation method, we denote the mean and variance of the solar azimuthal angle estimation utilizing the feature y by µ ap , σ ap . We denote the mean, variance and mean error of the solar elevation angle estimation utilizing features y by µ ep , σ ep , e ep . The calculated results of the four experiments are reported in Fig. 12 (a), where we find that variances σ ap and σ ep in all experiments have relatively low values (less than 3°), which corresponds with the practical conditions of experiment E i . Observing the mean error term e ep , we find that the mean error values of E 1 , E 2 , E 4 are acceptable, while the error value of E 3 is relative large 4 . We claim that such undesirable estimation error occurs in practice since the proposed estimation method purely relies on image feedback for obtaining transformation matrices, which can be sometimes inaccurate. Such error can reduced if more sensors, e.g. compass and IMU, are integrated in the system.

Besides the position-based feature vector y, there are other types of sensor-based geometric features which can be used for solar orientation estimation. For example, when the lens plane is set parallel to the ground plane, the orientation of the principal axis of the ellipse (denoted by "orientation-based" feature in Fig. 11 b) is coupled with the solar azimuthal angle, and the eccentricity (shape) of the ellipse is coupled with the solar elevation angle. The two geometric features for solar orientation estimation can be selected arbitrarily as long as they are not linearly dependent. Yet, we find that the orientation and shape features of the light spot are more sensitive to inaccurate system calibration compared to the position-based features. For example, we denote the variance of the solar azimuthal angle estimation using the orientationbased feature by σ ao . The comparison between the σ ao and σ ap is depicted in Fig. 12 (b). The position-based feature that yields more stable estimations is finally selected for the proposed method.

C. Robotic Solar Painting Experiments

We conducted four tests to evaluate the proposed system's accuracy of focusing solar energy to the given target positions. We denote the four tests by T i=1,2,3,4 , and the k th target point in T i is denoted by w p k i = [x k i , y k i ] . For each test, the corresponding histograms are depicted in Fig. 13 (c). In general, the accuracy performance of the developed system is sufficient for conducting solar painting tasks. We conducted solar painting experiments with various target patterns and different scaling factors to validate the proposed methods. The process and results of four representative experiments are demonstrated in Fig. 14. Depending on the amount of target positions sent to the robot's motion controller, the operation time for each heliography experiment varies from 7 minutes to 15 minutes. One complication of the proposed robotic approach comes from the difficulty to accurately detect the ground plane. In our method, we assume the ground is flat such that the homography transformation is valid and can be determined by an ArUco marker. However, as the wooden surface may (slightly) deform after being exposed to solar energy for a long period, this constraint may not be exactly satisfied. Also note that as many commercial depth cameras rely on infrared light, thus, they may perform poorly in outdoor environments. Another limitation of the proposed method is the difficulty to observe the point where solar energy is concentrated. One possible solution is to calibrate the thermal camera [START_REF] Zhang | Towards a multispectral rgb-ir-uv-d vision system-seeing the invisible in 3d[END_REF] that is invariant to light intensity for determining the location of this point.

The accompanying video demonstrates the performance of the proposed robotic heliography system with multiple experiments: https://github.com/JinliBot7/SR/raw/main/video.mp4

VI. CONCLUSION

In this paper, we presented a novel sensor-based robotic system that is capable of automatically performing heliography, an art form that involves the creation of pictorial patterns on a wooden canvas by concentrating solar energy over it. In order to focus the sun rays, the robot's end-effector rigidly grasps a Fresnel lens, and is equipped with thermal and RGB cameras. To perform the procedure, a desired input figure is analyzed to extract contour features, which result in a series of target points where solar energy will be concentrated. Then, the motion controller calculates the required robot position such that solar energy is concentrated over each of these points, The sun direction is estimated in real time based on visual feedback to maximize the sun rays passing through the lens.

To create carbonized points of uniform size, our sensor-based method monitors the wood pyrolysis and thus calibrate the focal distance to the surface As a preliminary step towards robotic painting by harnessing solar energy, our work has achieved promising results. Specifically, we identify that the key to the success of the experiments lies in the accurate estimation of the sun direction, as well as on proper determination of the lens-plane distance. To estimate the sun direction, we design an iterative update rule that is based on the interaction (Jacobian) matrix of the system, which captures the relationship between the sun angle increments and the feature variations. We select a positionbased feature vector over an orientation-based feature, since it exhibits smaller variance and thus leads to more accurate and stable estimation. The effectiveness of our innovative heliography system is verified by autonomously drawing various representative patterns, namely a university logo, letters, a human portrait, and numbers.

There are various limitations associated to our current approach. For example, our system cannot function well during cloudy or rainy days, although this issue is also encountered by other solar energy-dependent applications [START_REF] Li | Quantitative evaluation of the impact of cloud transmittance and cloud velocity on the accuracy of short-term dni forecasts[END_REF]. When the sunlight is completely blocked by clouds, the temperature will not rise to the critical point and the robot will wait at a target position until the sunlight reappears. When the sun is partially blocked, pyrolysis cannot be reached in a consistent manner, see e.g., Fig. 15(d). On the other hand, working for a long period under intensive sunlight could also cause problems, since overheating the motors and other electrical components might bring potential risks to the robot. Since adding a heat shield might hinder robot motion, the current solution is to manually place the main body of the system under the shadow and only expose to sunlight the end-effector part which carries the lens. Yet, the position of the shadow moves as time passes, so it may cover the end-effector after a period of time, causing the termination of an ongoing task. Besides, since our robotic prototype is a fixed-base manipulator, the limited workspace of the robot end-effector limits the size of the heliographic painting.

For future work, there are several directions along which we are interested in investigating. A feasible solution to paint a large picture is to use a mobile/wheeled a robot manipulator such that the painting range can be considerably extended [START_REF] Hu | Nonlinear model predictive control for mobile medical robot using neural optimization[END_REF]. For motion generation, it is worth investigating the transfer of drawing skills from a human teacher to a robot, so that the procedure of image processing the desired pattern can be bypassed [START_REF] Zahra | A bio-inspired mechanism for learning from mirrored demonstrations[END_REF]. Besides, proper coordination between the end-effector and joint links would be necessary when the robot moves into a configuration where the links block the sunlight [START_REF] Duan | Learning to sequence multiple tasks with competing constraints[END_REF]. For the sun angle estimation, we would like to integrate the recent advances in the determination of plane of array (POA) solar irradiation on inclined surfaces, so that the optimal tilt and azimuth angles could be obtained to maximize the power output [START_REF] Li | Spectral solar irradiance on inclined surfaces: A fast monte carlo approach[END_REF].

Fig. 1 .

 1 Fig. 1. (a) Illustration of heliography performed by an artist [2] and (b) by a robotic arm.

Fig. 2 .

 2 Fig. 2. Workflow of the proposed robotic painter.

Fig. 3 .Fig. 4 .

 34 Fig.3. A schematic of the system with relevant frames. The system is composed of a 3-DOF robot arm with two additional motors (controlled by a Raspberry Pi), manipulating a customized end-effector. The end-effector carries a spherical Fresnel lens, an RGB camera, and a thermal camera. An ArUco marker defines the World Frame.

  (a)). For an incident point o p l1 on the convex boundary, and an incident point o p l2 on the flat boundary, the corresponding active surface normals are (see Fig. 4(b)):

Fig. 5 .

 5 Fig. 5. Illustration of the geometrical characteristics of the projected points. Remember to change r l to r l i

Fig. 6 .

 6 Fig.6. The effect of the solar orientation variation on the light spot when the lens plane is parallel to the ground plane.

Fig. 7 .

 7 Fig. 7. Simulations of the process of lens reorientation and focus, with polar angle ϕ 0 = 5/6π and azimuthal angle θ 0 = 1/6π

Fig. 8 .

 8 Fig. 8. Experiments designed to determine the parameters for the motion control algorithm. (a) Solar concentration experiments with fixed duration and varying lens-plane distance. (b) Representative thermal images captured before, during and after the wood pyrolysis process; (c) Variation of the pyrolysis region with respect to time; (d) The size of a carbonized point can be controlled by the selection of Np.

Fig. 9 .

 9 Fig. 9. Feature extraction and scaling from target figures for path planning. The obtained target points were used for conducting real solar painting experiments as reported in Sec. V-C.

Fig. 10 .

 10 Fig. 10. Image-based focal length calibration experiment. (a) includes snapshots of captured light spot and extracted contour; (b) shows the corresponding lens-plane configuration in the developed simulation; (c) shows the comparison between the light spot area from the real images and the simulation.

Algorithm 1 :

 1 Motion control of the robot given P w Robot-World Calibration ← Section II-B; Initial Solar Orientation Estimation ϑ ← Section III-C; Target Position List ← P w = { w p * 1 , . . . , w p * k }; Target Position Index i ← 1; while i ≤ k do while N p ≤ N p do N p ← Thermal Camera; d lp = d lp ; r p * l ← (14); end while N p ≥ 0 do N p ← Thermal Camera; d lp ← d lp + δ d lp ; r p * l ← (14); end

Fig. 11 .

 11 Fig. 11. Experiments evaluating the performance of the solar orientation estimation utilizing different geometric features. (a 1 )-(a 4 ) include images captured from four experiments with varying lens yaw angles, and (c 1 ) -(c 4 ) show the corresponding robot configuration from the top view. (b) annotates the components in the robot configuration visualization.

Fig. 12 .

 12 Fig. 12. A quantitative evaluation of the performance of different estimation methods. The position-based method outperforms the orientation-based method, in terms of accuracy and stability.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Validation of the positional accuracy of the developed system. The mean of the horizontal error is 0.099 cm, and the mean of vertical error is 0.061 cm.

Fig. 15 .

 15 Fig. 15. Examples of failed solar paintings, which are caused by (a) inappropriate feature selection in the path planning stage, (b) the lens being shadowed as the solar orientation changed, (c) inaccurate solar orientation estimation and (d) a cloudy situation when Np is set too small.

The sun elevation angle is the angular height of the sun, which is equal to 0 at sunrise. The sun azimuthal angle is the angle between the projection of sun's centre onto the horizontal plane and the due north for the northern hemisphere[START_REF] Kalogirou | Solar thermal systems: Components and applicationsintroduction[END_REF]. The solar orientation under this conventional coordinate system can be easily obtained given the longitude, latitude and time information. However, as the orientation of the robot frame with respect to north is unknown, the solar orientation under the robot frame requires further estimation.

The feature vector y is actually defined in frame (xyz)o. As the transformation between (xyz)o and (xyz)r is known, y can be easily computed from the red ellipse in (xyz)r.

[START_REF] Calinon | A humanoid robot drawing human portraits[END_REF] The error between the estimated and real solar elevation angle can be obtained since the robot and the ground plane are all placed on the ground. Yet, the real value of the robot's orientation with respect to north is unknown since there is no external device integrated with the system. Thus, we use the variance as an alternative to evaluate the estimation of solar azimuthal angle.

This work is supported in part by the Research Grants Council of Hong Kong under grant 15212721 and in part by the Jiangsu Industrial Technology Research Institute Collaborative Funding Scheme under grant ZG9V.