Method for estimating an orientation change characteristic

Abstract

A method for estimating an orientation change characteristic of an object from images depicting the object captured at different times. For a plurality of candidate values of the orientation change characteristic, a first region in a first image and a second region in a second image are identified, wherein the first region and the second region would represent a same portion of a surface of the object if the orientation change characteristic was equal to the candidate value. For the plurality of candidate values, a value of a cost function is determined, based on intensity values of a pixel in the first region and a pixel in the second region, and a parameter set comprising a respective value for one or more illumination parameters. The parameter set represents illumination conditions of the object. The orientation change characteristic is estimated based on the values of the cost function.

Description

[0001]METHOD FOR ESTIMATING AN ORIENTATION CHANGE CHARACTERISTIC Technical Field The present invention relates to a method, processor, and non-transitory computer readable medium for estimating an orientation change characteristic of an object from a plurality of images depicting the object captured at respective times. In examples disclosed herein, example methods and processors are used to capture video of a sporting activity, such as a sporting event or a practice session. During a sports event, it can be desirable to know how fast a ball is spinning. Some methods for estimating the angular velocity of a ball from images depicting the ball have been developed. Such methods generally rely on a model of the appearance of the surface of the ball. By comparing the appearance of a ball in an image with the known appearance of the ball, the orientation of the ball at the time at which the image was taken can be determined. Then, the orientation of the ball calculated at two different times can be used to determine the angular velocity of the ball. As mentioned, such methods rely on knowledge of the appearance of the ball. This means that, for balls that look different to each other, bespoke models are required. Each bespoke model must include, for every point on the surface of the ball, a value indicating the colour of that point. Furthermore, in order to ensure that the camera reproduces an image of the ball that matches the model, the camera must be carefully colour-calibrated. Summary According to a first aspect of the present invention, there is provided a method for estimating an orientation change characteristic of an object from a plurality of images depicting the object captured at respective different times, the method comprising: for each of a plurality of candidate values of the orientation change characteristic for the object: identifying a first region in a first image of the plurality of images and a second region in a second image of the plurality of images, wherein the first region and the second region would represent a same portion of a surface of the object if the orientation change characteristic was equal to the candidate value; and determining, based on pixel intensity values of one or more pixels in the first region, pixel intensity values of one or more pixels in the second region, and an illumination parameter set comprising a respective value for each of one or more illumination parameters, the illumination parameter set representing illumination conditions of the object in the plurality of images, a value of a cost function; and estimating the orientation change characteristic based on the values of the cost function. Optionally, identifying the second region comprises identifying the second region based on the first region and a shape of the object. Optionally, determining the value of the cost function comprises: determining, based on the illumination parameter set, at least one illumination-adjusted difference between the pixel intensity values of the one or more pixels in the first region and the pixel intensity values of the one or more pixels in the second region; and determining the value of the cost function based on the illumination-adjusted difference. Optionally, the illumination parameter set is a first candidate illumination parameter set of a plurality of candidate illumination parameter sets, and the method comprises: for each of the plurality of candidate values of the orientation change characteristic for the object: for each of the plurality of candidate illumination parameter sets, determining, based on the pixel intensity values of the one or more pixels in the first region, the pixel intensity values of the one or more pixels in the second region, and the candidate illumination parameter set, an illumination-specific value of the cost function; and estimating the orientation change characteristic based on the illumination-specific values of the cost function. Optionally, the method comprises determining the illumination parameter set representing the illumination conditions of the object in the plurality of images by selecting a candidate illumination parameter set, from among the plurality of candidate illumination parameter sets, based on the values of the cost function. Optionally, estimating the orientation change characteristic comprises selecting a candidate value of the orientation change characteristic, from among the plurality of candidate values, that has a lowest value of the cost function among the values of the cost function. Optionally, the illumination parameters comprise at least one of: a location of a light source relative to the object, a direction of illumination of the object, an intensity of the illumination, a colour of the illumination, a specular reflection property of the object, a diffuse reflection property of the object, or an ambient reflection property of the object. Optionally, the orientation change characteristic is a first angular velocity for a first time, and estimating the first angular velocity comprises, for each of the plurality of candidate values of the first angular velocity: determining the value of the cost function further based on a difference between the candidate value of the first angular velocity and a second angular velocity of the object for a second, different, time, so as to reduce a difference between the estimated first angular velocity and the second angular velocity. Optionally, determining the value of the cost function comprises weighting the difference between the candidate value of the first angular velocity and the second angular velocity according to a difference between the first time and the second time. Optionally, identifying the second region comprises: determining a first position of the portion of the surface of the object in 3D space, from the first image; determining, based on the first position and the candidate value of the orientation change characteristic, a second position of the portion of the surface of the object in the 3D space, the second position representing a position of the portion of the surface of the object in the 3D space at a time of capture of the second image; and determining, based on the second position of the portion of the surface of the object in the 3D space, the second region in the second image. Optionally, determining the first position comprises determining the first position based on a known size of the object and based on a measured size of the object in the first image, and / or identifying the second position comprises identifying the second position based on the known size of the object and based on a measured size of the object in the second image. Optionally, the orientation change characteristic comprises an angular speed, at least one component of angular velocity, or a direction of angular velocity, of the object. Optionally, the object is a sports ball. Optionally, the method comprises: determining, based on a time difference between the first image and the second image, and on a maximum expected angular speed of the object, that a maximum expected angular displacement of the object between the first image and the second image is lower than a threshold, wherein identifying the first region in the first image and the second region in the second image is performed in response to determining that the maximum expected angular displacement is lower than the threshold. According to a second aspect of the present invention, there is provided a processing system configured to perform the method of the first aspect. According to a third aspect of the present invention, there is provided a system comprising: the processing system of the second aspect, and one or more cameras, the one or more cameras being configured to provide the plurality of images to the processing system. According to a fourth aspect of the present invention, there is provided a non- transitory computer readable medium comprising instructions which, when executed by a processor, cause the processor to perform the method of the first aspect. Figure 1 shows a schematic diagram of a method according to the invention; Figure 2 shows a schematic diagram of a system for implementing the method; Figure 3 shows a perspective view of the system of Figure 2 when deployed at a basketball game; Figure 4 shows a schematic geometric diagram indicating the relative arrangement of the components of the system when used to implement the method; Figure 5 shows schematically a technique for determining the position of an object using two cameras; Figure 6 shows a schematic illustration of a technique for determining corresponding regions in different images, as used in the method; and Figure 7 shows an example technique for estimating an orientation change characteristic, according to the method. Detailed Referring to Figure 1, there is illustrated a method 100 for estimating an orientation change characteristic of an object 50 from a plurality of images 40, 41 depicting the object 50 captured at respective times. In some embodiments, the orientation change characteristic is an angular velocity of a sports ball. The sports ball is in play in a sports match and spins as it moves through the air. Images of the ball are captured at different times by a rig of cameras 110(a)-(g) positioned around a playing space, for example a pitch or court on which the sports match is being played. In broad overview, the method 100 comprises: - in step 102, for each of a plurality of candidate values of the orientation change characteristic for the object 50, identifying a first region x0 in a first image 40 of the plurality of images 40, 41 and a second region x1 in a second image 41 of the plurality of images 40, 41, wherein the first region x0and the second region x1 would represent a same portion of a surface of the object 50 if the orientation change characteristic was equal to the candidate value; - in step 104, for each of the plurality of candidate values of the orientation change characteristic for the object 50, determining, based on pixel intensity values of one or more pixels in the first region x0, pixel intensity values of one or more pixels in the second region x1, and an illumination parameter set comprising a respective value for each of one or more illumination parameters λ1, λ2, …, λk, the illumination parameter set representing illumination conditions of the object 50 in the plurality of images 40, 41, a value of a cost function; and - in step 106, estimating the orientation change characteristic based on the values of the cost function. Accordingly, a method 100 for estimating an orientation change characteristic of an object 50 is provided for. Broadly speaking, the method 100 works by determining the candidate value of the orientation change characteristic (e.g. angular velocity of the object 50) that is consistent with the observed changing appearance of the object 50 across the plurality of images 40, 41. As an introductory illustrative example, assume that the surface of the object 50 is painted in greyscale according to a pattern in which no part of the surface has the same greyscale tone as any other part of the surface, the object 50 is uniformly illuminated by planar radiation from all angles at all times, and the object 50 is a perfect diffuse reflector of light. It is “guessed” that the object 50 rotates through a particular angle, e.g.3 degrees, between the images. A pixel intensity value in the first image 40 representing a portion of the surface of the object 50 is chosen. The object 50 is “rotated” by 3 degrees to obtain the theoretical (3D) location at which this portion of the surface would be in the second image 41. If the pixel intensity value in the second image 41 representing this theoretical location is similar to the pixel intensity value in the first image 40, then it is likely that the object 50 did indeed rotate by approximately 3 degrees between the images. Meanwhile, if the pixel intensity values are very different, then it is likely that the pixel from the first image 40 represents a different part of the surface of the object 50 to the pixel in the second image 41, and hence the assumption that the object 50 rotated by 3 degrees was incorrect. This information regarding the accuracy of the assumed angular displacement can be reflected in the cost function, by setting the cost function equal to the difference between the pixel intensity values. If the value of the cost function is low, then it is likely that the assumed angular displacement is correct, while if it is high, then it is likely that the assumed angular displacement is incorrect. The cost function can be computed for various different “guessed” rotation angles. To estimate the angular displacement, the rotation angle corresponding to the lowest value of the photogrammetric cost function can be identified. By determining a cost function for each of a plurality of candidate values of the orientation change characteristic (e.g. angular velocity, or angular displacement), the orientation change characteristic can be estimated without knowledge of the appearance of the surface of the ball. That is, instead of identifying the orientation of the ball at each time by looking for known markings on the ball in each image, the orientation change characteristic can be estimated, irrespective of the appearance of the surface of the ball, by determining the angular velocity that most closely matches the observed pixel intensity values. In other words, the appearance of part of the surface of the ball (or a portion thereof) facing the camera 110 can be gleaned from the first image 40, and the same can be done for the second image 41; by comparing the appearance of a portion of the surface of the ball that would be present in both images assuming a particular angular velocity, it can be determined whether the assumed angular velocity is correct. The above example assumes that the object 50 is uniformly illuminated from all angles at all times. Of course, in reality, objects are often illuminated in a non-uniform fashion, and reflect light in a non-uniform fashion. Thus, in the above example, even if the assumed angular velocity matches the true angular velocity, the pixel intensity values in different images representing the same portion of the surface of the object 50 would in general not match each other, because the object 50 is rotating and thus this portion of the surface of the object 50 may be more or less illuminated in the second image 41 than the first image 40. To address this, the method 100 comprises determining values of the cost function using not only a plurality of candidate values of the orientation change characteristic, but also an illumination parameter set representing the illumination conditions of the object 50 in the plurality of images 40, 41. This allows the orientation change characteristic to be estimated reliably and accurately in a range of different lighting conditions and for a range of objects with different reflection properties. System overview Referring to Figure 2, there is illustrated a system 200 according to an example. The system 200 comprises a server 210 (also referred to as a processing system) and a plurality of cameras 110(a)-(g). The method 100 comprises receiving the plurality of images 40, 41, at the server 210, from the plurality of cameras 110(a)-(g). The cameras 110(a)-(g) are positioned at respective locations in or near a sporting event. The server 210 comprises an input interface 213, an output interface 214, a processor 211, and a memory 212. The processor 211 and the memory 212 are configured to perform the method. The memory 212 stores instructions which, when executed by the processor 211 cause the processor 211 to perform the method. The instructions may be stored on any computer readable medium, for example any non- transitory computer readable medium. Referring to Figure 3, there is illustrated an example of the system 200 of Figure 2 as deployed to analyze a basketball game. Each of the cameras 110(a)-(g) shown in Figure 3 has a corresponding, different viewpoint 114 from which it captures video. As shown, most of the cameras 110(a)-(g) are arranged at various fixed positions and orientations around the basketball court. In some examples, some of the cameras 110(a)-(g) could be held by coaches or fans (e.g., the fans’ own personal devices could function as the cameras 110(a)-(g) of the system 200). Furthermore, although a basketball game is depicted in Figure 3, it will be understood that this is merely illustrative and that the system 200 of Figure 2 is suitable for deployment at many other types of sporting activity and, indeed, is suitable for deployment in non-sporting environments, such as a non-sporting live event (e.g., a concert, a comedy show, or a play) or a non-sporting practice session (e.g., a music practice, or a rehearsal for a play). Returning now to the system 200 of Figure 2, it should be noted that, in some examples, the server 210 may be located in the same physical location as the portable electronic devices. For example, in a situation where the system 200 is deployed at a sporting activity, for instance as illustrated in Figure 3, the server 210 may be located in a server room at the venue where the sporting activity is taking place. Alternatively, the server 210 might be located in a truck parked on-site at the venue. However, in still other examples, the server 210 could be a remote / cloud server. According to the method 100, the cameras 110(a)-(g) have been synchronised to capture frames at the same time, for example according to the techniques described in the patent application with publication number PCT / US2023 / 036480. Alternatively, the server 210 may instruct the cameras 110(a)-(g) to capture images at predefined, future, timestamps. The cameras 110(a)-(g) continuously capture frames depicting a sports ball that is in play in the sporting activity. The sports ball may be, for example, a basketball, a volleyball, a tennis ball, a soccer ball, or a cricket ball. Each image is an RGB image of resolution H x W, associated with a timestamp t stored by the respective camera 110 (a)-(g) and transmitted to the server 210 together with the image. Ball detection The method 100 comprises, for each of the plurality of images 40, 41, detecting a portion of the image in which the sports ball is depicted. Various techniques for detecting a ball in an image are known in the art, but for completeness, we describe techniques for ball detection that may be used in implementing the methods disclosed herein. A neural network, trained to detect an object with a particular shape (e.g. a circular object representing a basketball, or a non-circular object such as a rugby ball) in an image, is used to obtain an initial estimate of the portion of the image in which the ball is depicted, as well as the diameter of the circle (if the ball is circular). Subsequently, a square region (if the ball is circular) of the image having side length greater than the diameter of the circle (e.g. 1.2 times the diameter of the circle) is cropped, and a Hough transform is applied to determine a refined estimate of the portion of the image in which the ball is depicted. The portion of the image in which the ball is depicted is referred to as a mask, and comprises the pixels that represent the ball. Alternatively, instead of using the neural network to obtain the initial estimate, a Hough transform may be applied directly to the entire image to obtain the position of the ball. 3D ball position calculation Once the ball has been identified in a given image, the 3D position of the ball in space can be determined, as illustrated with reference to Figure 4. Again, various techniques for this exist in the art, but for completeness, we describe techniques for determining the 3D position of the ball that may be used in implementing the methods disclosed herein. Firstly, a camera calibration process can be performed to determine extrinsic and intrinsic parameters of the camera 110. Extrinsic parameters include the position (Xc, Yc, Zc) (in a 3D space, for example the position (Xc, Yc, Zc) of the camera 110 relative to a pitch on which the sporting activity is played) and orientation of the camera 110. The orientation may be represented using a quaternion. Alternatively, the orientation may be represented using a vector whose direction is parallel to the ray direction corresponding to the centre of the field of view of the camera 110. Alternatively, the orientation may be represented using a rotation matrix, where the first column of the matrix is a vector whose direction is parallel to a vertical direction y (see e.g. Figure 4) in the plane of an image 40, 41 captured by the camera 110, the second column of the matrix is a vector whose direction is parallel to a horizontal direction x in the plane of the image 40, 41 captured by the camera 110, and the third column is parallel to the ray direction corresponding to the centre of the field of view of the camera 110. Here, the vertical direction y and the horizontal direction x in the image 40, 41 are directions that are mutually perpendicular and perpendicular to the ray direction corresponding to the centre of the field of view of the camera. A camera calibration process for determining the extrinsic parameters is described in patent publication US 10600210 B1. Alternatively, the position (Xc, Yc, Zc) of the camera 110 may be determined using a GPS sensor in the camera 110 (not shown). Intrinsic parameters include the focal length f of the camera 110, the field of view of the camera 110, and parameters describing optical distortions associated with the camera 110. The intrinsic parameters may also be determined according to the camera calibration process described in patent publication US 10600210 B1; or alternatively may be obtained from a memory of the camera 110 or the memory 212 of the server 210. The intrinsic parameters, the diameter Dmof the ball measured in the image, and the known diameter Dr of the ball, can be used to determine the distance L between the centre of mass of the ball and the camera 110, by inverting the following equation: , is the distortion ratio of a distortion function 2D pixel location x (where the pixel location x = (0, 0) represents the centre of the image). The distortion ratio is equal to 1 in the case of no distortion. The extrinsic parameters, the intrinsic parameters, the position (xb, yb) of the centre of the ball measured in the image (i.e. the pixel representing the centre of the hemispherical surface of the ball facing the camera 110), and the distance L between the ball and the camera 110, can then be used to determine the 3D position (Xb, Yb, Zb) of the centre of mass of the ball. Specifically, the position (xb, yb) of the centre of the ball measured in the image can be used together with the known intrinsic parameters to determine the unit vector ^^for the line adjoining the camera 110 and the centre of mass of the ball. The 3D position (Xb, Yb, Zb) of the centre of mass of the ball can then be calculated according to equation (1).(^^, ^^, ^^)=(^^, ^^, ^^)+ ^ × ^^(1)The 3D position (Xb, Yb, Zb) of the centre of mass of the ball may be calculated from multiple images taken from different cameras at the same time, to improve accuracy. For instance, for each of a plurality of candidate 3D positions, a reprojection error can be calculated for each image. The reprojection error for a given 3D position and image is, for example, the angle between the line ^^(calculated for that image using the method described above) adjoining the corresponding camera to the centre of mass of the ball, and the line adjoining the camera to the candidate 3D position. The reprojection error for a given 3D position and image may instead be the 2D distance (measured e.g. in pixels) between the point in the image representing the centre of mass of the ball (calculated using the method described above) and the point in the image representing the candidate 3D position. Alternatively, the reprojection error for a given 3D position and image may be the 3D distance (measured e.g. in metres) between the 3D position of the centre of mass of the ball (calculated using the method described above) and the candidate 3D position. For a given candidate 3D position, the reprojection error is calculated for each image and summed to obtain a total reprojection error. The candidate 3D position with the lowest total reprojection error is selected as the 3D position of the centre of mass of the ball. The points on the surface of the ball satisfy(X − X^)^+(Y − Y^)^+(Z − Z^)^=(^^ / 2)^. (For non-spherical objects, an appropriate equation describing the known shape of the object can be used in place of this equation; thus, more generally, identifying the second region x1comprises identifying the second region x1based on the first region x0and a shape (e.g. a known shape) of the object 50. Further details regarding such equations are provided below.) For a given such point(^, ^, ^), the line adjoining this point to the camera 110 can be determined as(^, ^, ^)−(^^, ^^, ^^). This line can be used together with the known intrinsic parameters to determine the pixel (x, y) representing this point on the surface of the ball in the image. By extension, for a particular point on the surface of the ball, the pixel (x, y) representing this point on the surface of the ball in an image taken by any camera can be determined using the extrinsic and intrinsic parameters of the camera. The conversion of a point on the surface of the ball X = (X, Y, Z) to a pixel x = (x, y) is denoted as x = ϕc(X), where c denotes the camera used to capture the image containing pixel x. The converse of this, in which the point on the surface of the ball X = (X, Y, Z) represented by a pixel x = (x, y) is found, involves identifying the intersection of the line ^^for the pixel x with the surface of the ball; the line is denoted as In general there will be two intersection points; the intersection point closer to the camera 110 is retained and the result is denoted ζ(ϕ -1 c (x)). Such conversion methods may be applied using a variety of camera models, such as a pinhole model or an RPC model. The above describes a method for calculating the 3D position (Xb, Yb, Zb) of the centre of mass of the ball from one image. However, if two cameras with known extrinsic and intrinsic parameters capture an image of the ball at the same time, the 3D position (Xb, Yb, Zb) of the centre of mass of the ball may be determined by simple triangulation, without knowledge of the size of the ball. That is, the intersection point of the line ^^1adjoining a first camera to the centre of mass of the ball, with the line ^^2adjoining a second camera to the centre of mass of the ball, can be found, and used as the 3D position (Xb, Yb, Zb) of the centre of mass of the ball, as illustrated in Figure 5. As mentioned, the object 50 does not need to be spherical; the ball may be for example a rugby ball. More generally, the geometry G of the object 50 can be represented by a solid geometric object modelled by any of the following: 1. An implicit equation 2. An explicit 3. A mesh wherein is the set of triangles representing the mesh, with an intersection function which can be used to determine the intersection point between a 3D line (such as ^^1) and the geometry G. Such geometric models are per se well known in the art and will not be discussed further herein. Furthermore, the determination of the 3D position (Xb, Yb, Zb) of the centre of mass of the ball is not limited to the above methods; for example, the position may be received from a GPS sensor embedded within the ball. Image selection Once the 3D position (Xb, Yb, Zb) of the ball (or more generally object 50) has been calculated, the method 100 comprises, for at least one of the times at which images are captured, selecting a subset of the images captured at that time, to be used to estimate the orientation change characteristic (e.g. angular velocity), based on distances between the respective cameras 110(a)-(g) and the ball (or the size of the ball as shown on the image), and / or on the resolutions of the cameras 110(a)-(g). The cameras 110(a)- (g) may be numerous (for example, there may be up to 84 cameras), and for one given timestamp, a given part of the surface of the ball may be depicted by images from several different cameras 110(a)-(g). To reduce the processing requirements, a subset of cameras is chosen appropriately. Firstly, for each point (e.g. a small surface area) on the surface of the ball, and each camera 110 whose image includes that area, the number of pixels depicting that area is calculated, based on the distance L between the ball and the camera 110. If the image from more than one camera depicts that point, then the camera whose image depicts the highest density of pixels at that point is chosen. This reduces the number of representations of (i.e. pixels representing) each point on the surface of the ball, to one. In general, because some of the cameras are closer to the ball than others at a given capture time, the number of cameras whose images are used to represent the ball for a given capture time will be substantially lower than the total number of cameras. In any case, it is generally sufficient to use, for a given capture time, four images depicting the ball, and hence in method 100, a maximum of four images are used. It should be noted that the image selection technique described above is merely optional; in principle, all that is needed to estimate the orientation change characteristic is two images depicting the object 50 captured at different times. Image pairing Once the images for each capture time have been selected, images taken at different times can be paired in order to determine the cost function. Broadly speaking, the method 100 comprises: determining, based on a time difference between the first image 40 and the second image 41, and on a maximum expected angular speed of the object 50, that a maximum expected angular displacement of the object 50 between the first image 40 and the second image 41 is lower than a threshold, wherein identifying the first region x0 in the first image 40 and the second region x1 in the second image 41 is performed in response to determining that the maximum expected angular displacement is lower than the threshold. Consider a pair of images depicting the object 50, taken from the same camera 110 at times t0 and t0 + δt. In general, the object 50 will have rotated in the interval δt between the two images. By using a maximum expected angular speed wmax, the maximum expected angular displacement can be calculated, as wmax*δt. The maximum expected angular speed may be obtained empirically. For example, where the object 50 is a sports ball, such as a basketball, a maximum recorded angular speed of a basketball can be used as the maximum expected angular velocity. The maximum recorded angular speed may have been estimated using other methods, or using the methods described herein with a high-framerate camera. By ensuring that δt is not too large, it can be ensured that the ball has not rotated by more than half a revolution during the interval δt. Specifically, the maximum permitted time interval between two images to be compared is δtmax = π / wmax (with wmax measured in radians per second). (Nevertheless, in embodiments where the orientation change characteristic is not the angular velocity but instead the difference in orientation of the ball between the times at which two images are captured, it is not necessary to impose this requirement.) In practice, it is desirable to ensure that a sufficiently large part (i.e. a threshold percentage, such as 15%, or alternatively any value in the range 1-40%, or 10-20%) of the surface of the ball appears in both images, to improve the quality of the comparison between the images (the comparison between images is described below under “Cost function”). Therefore, it is desirable to ensure that the ball has rotated by a maximum of a predefined angle σ which is less than π radians. This threshold reduces the permitted time interval between two images to a value ε. From a given camera, any two images taken by that camera with a time separation less than ε can constitute a pair of comparable images. In a multi-camera setup, it is possible to compare images that are taken by different cameras (e.g. cameras with different orientations). In such a case, the permitted time interval between the images is lower than ε, because it must be ensured that the ball has rotated by less of an angle (in order to ensure that a sufficiently large part of the surface of the ball appears in both images). In the worst case scenario, the ball rotates with angular speed wmax in a direction opposite to θc2 - θc1, where θc1 is the orientation of the camera C1 used to capture the earlier image in the pair, and θc2is the orientation of the camera C2 used to capture the later image in the pair. In such a case, the maximum permitted time δtmax satisfies wmax*δtmax + |θc2 - θc1| = σ. Since the rotation axis of the ball is in general unknown, the maximum permitted time between a compared pair of images taken by cameras C1 and C2 is (σ - |θc2- θc1|) / wmax. Pairs of images can thus be selected appropriately. As mentioned, for each time t, four images are selected depicting the ball. Therefore, for a pair of times t0, t1, up to four pairs of images are selected. Where the cameras are synchronised, a pair of images may comprise a pair of adjacent frames. In principle, however, only one pair of images, captured at different times, is required. The method 100 does not rely on the above techniques for 3D ball position calculation, image selection, and image pairing; for example, where the same camera 110 is used to capture the first image 40 and the second image 41, the method 100 may begin at region pairing, which will now be described. Once appropriate image pairs have been selected, a first image 40 and second image 41 within a given pair can be compared. Referring again to Figure 1, as mentioned, the method 100 comprises, at step 102, for each of a plurality of candidate values of the orientation change characteristic (e.g. angular velocity) for the object 50 (e.g. a sports ball), identifying a first region x0in a first image 40 of the plurality of images 40, 41 and a second region x1 in a second image 41 of the plurality of images 40, 41, wherein the first region x0 and the second region x1would represent a same portion of a surface of the object 50 (e.g. sports ball) if the orientation change characteristic (e.g. angular velocity) was equal to the candidate value. As will now be described, identifying the second region x1comprises: determining a first position of the portion of the surface of the object 50 in the 3D space, from the first image 40; determining, based on the first position and the candidate value of the orientation change characteristic, a second position of the portion of the surface of the object 50 in the 3D space, the second position representing a position of the portion of the surface of the object 50 in the 3D space at a time of capture of the second image 41; and determining, based on the second position of the portion of the surface of the object 50 in the 3D space, the second region x1in the second image 41. An example region pairing method 102a for identifying the first region x0 and the second region x1, assuming a candidate value wt0 of the orientation change characteristic, is shown in Figure 6. Figure 6 shows an example in which the same camera 110 is used to capture both the first image 40 (at time t0) and the second image 41 (at time t1), but as mentioned above, in other examples, different cameras are used for a pair of images. The angular velocity wt0is expressed in the same 3D basis (i.e. relative to the 3D space) as the orientation of the cameras 110(a)-(g). Firstly, the first region x0 is selected from the first image 40. (As mentioned later, the following process is repeated for all regions x0 in the first image 40 that depict the object 50.) The first region x0includes one pixel. The disclosure herein discusses examples in which regions comprising individual pixels are used, but it should be understood that a pair of corresponding regions (i.e. a region in one image and a region in another image that would represent the same portion of the surface of the object 50 if the orientation change characteristic was equal to the candidate value) can be made up of multiple pixels. The region pairing method 102a comprises identifying a first surface region (first position) X0in 3D space, the first surface region X0being the part of the surface of the ball represented by the first region x0 in the first image 40. The techniques described above with reference to Figure 4 can be used to identify the first surface region X0. The region pairing method 102a comprises determining a second surface region (second position) X1 in 3D space, where the second surface region X1 is the same part of the surface of the ball as X0but at the time t1at which the second image 41 is captured, assuming the candidate value wt0 of the orientation change characteristic. That is, the angular velocity is assumed to take a particular value wt0. Here, the angular velocity is assumed to have a constant value wt0during the interval t0to t1. As mentioned above, this is a small time interval (during which the object 50 rotates at most by half a revolution) and hence the assumption that angular velocity is unlikely to change substantially during this interval. Determining the second surface region X1 in 3D space is performed based on the first surface region X0 in 3D space and the candidate value of the orientation change characteristic. The matrix that rotates a 3D vector by angle wt0(t1- t0) around the axis wt0 is denoted as R(wt0(t1 - t0)) in Figure 6. Suitable matrices for performing such rotation transformations are well known in the art and hence will not be discussed in greater detail. It should be noted that, while the methods described herein refer to estimating the angular velocity wt0, other orientation change characteristics, such as the angular displacement Θ (equal to wt0(t1- t0)) of the object between t0and t1, can be estimated instead. Where the orientation change characteristic is the angular displacement Θ, it is not necessary to know the times t0 or t1: instead, the rotation matrix can be simply written as R(Θ). Determining the second surface region X1in 3D space additionally comprises (if the centre of mass of the ball is moving), as shown in Figure 6, adding a vector (T1– T0), where T1 represents the 3D position of the centre of mass of the ball at time t1, and T0represents the 3D position of the centre of mass of the ball at time t0. T1and T0can be determined respectively from images captured at times t1and t0by the plurality of cameras 110(a)-(g), as described above with reference to Figure 4. The region pairing method 102a comprises determining the second region x1 based on the second surface region X1in 3D space. This is performed using the intrinsic and extrinsic parameters of the camera 110, as described above with reference to Figure 4. In examples in which the camera 110 used to capture the first image 40 is different to the camera used to capture the second image 41, these would be the intrinsic and extrinsic parameters of the camera used to capture the second image 41. In Figure 6, the transformation of the 3D surface region X1 to the pixel location x1 is denoted as x1 = ϕc(X1). Thus, the second region x1 may be determined from the first region x0 and the candidate value of the orientation change characteristic. For objects that are not spherical, the geometry of the ball Gtat time t can be written more generally as (referring to the basic definitions of the geometry G provided under “3D ball position calculation”): 1. For an implicit equation, , and 2. For an explicit equation, , where XtB represents the position of the centre of mass of the ball at time t, and represents the orientation of the ball at time t. In some cases (e.g. under certain candidate values of the angular velocity), the second surface region X1 will be on the side of the object 50 that faces away from the camera used to capture the second image 41. In such a case, the region x0is disregarded for the purposes of calculating the cost function below. Cost function As mentioned, the method 100 comprises, at step 104, for each of the plurality of candidate values of the orientation change characteristic, determining, based on pixel intensity values of one or more pixels in the first region x0, pixel intensity values of one or more pixels in the second region x1, and an illumination parameter set comprising a respective value for each of one or more illumination parameters λ1, λ2, …, λk, the illumination parameter set representing illumination conditions of the object 50 in the plurality of images 40, 41, a value of a cost function. The cost function P includes a photogrammetric term, which is shown in equation (2) below. The photogrammetric term is the sum, over regions x0 in the first image 40, of a difference function Δ[a, b] of the intensity values I (x0) of the one or more pixels in the first region x0 and the intensity values I (x1) of the one or more pixels in the second region x1. (As described below with reference to equation (4), this may be generalised to include respective sums from a plurality of pairs of images.) The difference function Δ[a, b] (described in further detail below) is parameterised by the values for the one or more illumination parameters λ1, λ2, …, λk. The illumination parameters λ1, λ2, …, λkare described in further detail below but may include: a location of a light source relative to the object 50, an intensity of the illumination, a colour of the illumination, a specular reflection property of the object 50, a diffuse reflection property of the object 50, or an ambient reflection property of the object 50. wt0 is (as above) the candidate value of the angular velocity of the object 50 (e.g. sports ball) at time t0. As described above, the second region x1is determined using the angular velocity wt0, hence the dependence of the photogrammetric term on the angular velocity wt0. Also as described above, a particular first region x0 is ignored (i.e. makes no contribution to the sum) where the corresponding second surface region X1 is on a side of the ball that is not depicted in the second image 41. An “illumination parameter set”, as referred to herein, comprises a respective value for each of one or more illumination parameters λ1, λ2, …, λk. Different “candidate illumination parameter sets” have different sets of values for the one more illumination parameters. For example, a first candidate illumination parameter set may have a first set of respective values for the location of a light source relative to the object 50, the intensity of illumination etc., and a second candidate illumination parameter set may have a second set of respective values for the location of a light source relative to the object 50, the intensity of illumination etc., the second set differing from the first set in that at least one of the values is different. The cost function is evaluated for each of a plurality of candidate illumination parameter sets (and for each of the plurality of candidate values of the angular velocity wt0). That is, the locations of the light sources relative to the object 50, the direction of illumination of the object 50, the intensity of the illumination etc. (also written as λ1, λ2, …, λk) and the candidate angular velocity wt0, are each varied, and for each combination of candidate values, the photogrammetric term is computed. Figure 7 illustrates an example in which only the x-component (wt0)xof the angular velocity and the x-position λ1of the light source (relative to a fixed point in 3D space) are varied and all other variables are assumed to be known; for each combination of a candidate value for the x-component (wt0)x of the angular velocity with a candidate value for the x-position λ1 of the light source, the value of the cost function is computed, and shown in the corresponding cell of the table in Figure 7. More generally, however, the angular velocity wt0has three components (all independent variables in equation (2) that are to be determined), and there are multiple illumination parameters as λ1, λ2, …, λk. By considering a variety of possible illumination conditions, the orientation change characteristic (e.g. angular velocity) can be estimated accurately even when the actual lighting conditions and the appearance of the surface of the object 50 are unknown. However, in other examples, the values of the illumination parameters λ1, λ2, …, λk are known (e.g. from a prior estimation of the angular velocity and illumination parameters), and only the angular velocity wt0 is varied. Difference function We now describe the difference function Δ, firstly using an example in which each region x0 in the first image 40 includes just one pixel. In such an example, the difference function Δ is, to put it simply, the lighting-adjusted difference between the intensity value I (x1) of the pixel in the second region x1and the intensity value I (x0) of the pixel in the first region x0. As a first example, the difference function Δ may be the square, or the modulus, of the difference ^^(^^) − ^^(^^) between a lighting-independent intensity value ^^(^^) of the pixel in the second region x1and a lighting-independent intensity value ^^(^^) of the pixel in the first region x0. The lighting-independent intensity values ^^(^^), ^^(^^) are determined by mapping (as described in further detail below) the intensity values measured by the cameras I (x0), I (x1) to respective lighting-independent intensity values, using the values of the illumination parameters. For example, the lighting-independent intensity values may represent an intrinsic value representing the (RGB-summed) colour of the ball. As a second example, the difference function Δ may be the square, or the modulus, of the difference between the intensity value of a first of the two pixels, e.g. I (x1), and a lighting-adjusted intensity value of the other pixel. In this case, the lighting- adjusted intensity value is determined by mapping the intensity value I (x0) to an intensity value ^^(^^) that would be measured if the illumination condition was the same as the illumination condition under which the intensity value of the first of the two pixels I (x1) was measured. The above examples may be generalised to the case in which each region x0in the first image 40 includes more than one pixel, by taking the difference between the sum of the (appropriately lighting-adjusted) intensity values in the first region x0and the sum of the (appropriately lighting-adjusted) intensity values in the second region x1. In any case, the “intensity value for a pixel” may refer to the sum of the red, green and blue component intensity values for the pixel, or to another function of these values, such as a greyscale intensity value obtained from these components, or a ratio between two of the three components. Lighting model We now describe methods for calculating the lighting-independent and lighting- adjusted intensity values using known lighting models. In one example, the Phong reflection model is used to model the illumination of the object 50. It is assumed that the parameters of the Phong reflection model are constant throughout the capture of the plurality of images 40, 41. The illumination parameters λ1, λ2, …, λk are: the location of each light source m relative to the object 50 (and hence the direction of illumination of the object 50 by each light source m), intensity of the specular component is of each light source, intensity of the diffuse component id of each light source, ambient lighting ia, specular reflection constant ks of the object 50, diffuse reflection constant kd of the object 50, ambient reflection constant kaof the object 50, and shininess constant α of the object 50. The Phong reflection model provides equation (3) for computing the illumination Ip of a point on the surface of the object 50: Here, Lmrepresents the direction vector from the point on the surface towards a light source m (i.e. the illumination direction for that light source), N represents the direction vector of the normal to the surface, V represents the direction vector from the point on the surface to the camera 110, and Rmis equal to . (For the first region x0, i.e. the first surface region X0, N and V are known, and Lm is parameterised by the location of the light source m. For the second region x1, i.e. the second surface region X1, N and V are parameterised by the angular velocity wt0, and Lmis parameterised by both the location of the light source m and the angular velocity wt0.) For a given set of values of the illumination parameters (location of each light source, is, id, etc.), and a given angular velocity wt0, the illumination Ipcoming from regions X0and / or X1can be computed, as Ip(x0) and Ip(x1) respectively. Then, modelling the colour of the ball as the sum of its intrinsic colour (represented by the lighting-independent intensity value) and the contribution of the diffuse and specular components of the illumination (represented by Ip), the lighting-independent intensity values can be calculated as I(x0) - Ip(x0) and I(x1) - Ip(x1), so that the difference function is or . In the case that a lighting-adjusted intensity value is used instead, the difference function is identical; the lighting-adjusted intensity value for e.g. pixel x0 can be calculated as I (x0) - Ip (x0) + Ip (x1), this value being subtracted from the (non-adjusted) intensity value I (x1) for pixel x1 in the difference function. Other lighting models may be used, such as a Gouraud model. While, according to the above Phong model, the direction of illumination Lm is calculated from the relative position of the light source with respect to the object 50 and is hence not an independent parameter, in other models, it may be an independent parameter. Estimation of orientation change characteristic (e.g. angular velocity) As mentioned, the method 100 comprises, at step 106, estimating the orientation change characteristic (e.g. angular velocity) based on the values of the cost function. Estimating the orientation change characteristic comprises selecting a candidate value of the orientation change characteristic, from among the plurality of candidate values, that has a lowest value of the cost function among the values of the cost function. Additionally, the illumination parameter set that correctly represents the illumination of the object 50 is determined, from among the plurality of candidate illumination parameter sets, based on the values of the cost function. Specifically, the set (candidate value of X-angular velocity, candidate value of Y-angular velocity, candidate value of Z-angular velocity, candidate value of λ1, candidate value of λ2,…, candidate value of λk) that has the lowest value of the cost function, is selected. Irrespective of whether the selected candidate value of the orientation change characteristic, and / or the selected candidate illumination parameter set, has the lowest value of the cost function among the values of the cost function, the orientation change characteristic can be estimated, and optionally the illumination condition determined, using an iterative solving method such as the Levenberg-Marquardt algorithm, or alternatively, gradient descent. For the iterative solving method, an initial candidate angular velocity, such as (0, 0, 0), and an initial illumination parameter set (e.g. where each component is a random number) can be used as a starting point. The iterative method uses values for each component of the angular velocity capped between -wmax and +wmax, i.e. no value outside this range is considered. Various iterative methods for finding a local or global minimum value of a function of multiple variables are well known in the art, and these methods will not be discussed herein. In the example shown in Figure 7, it can be determined that the combination of x-component of the angular velocity (wt0)xand x-position of light source λ1which leads to the lowest value of the cost function, is (20 rad s-1, 2 m). Smoothing term In a simple example, the cost function includes only the photogrammetric term referred to above, and the angular velocity is calculated for just one time. In some cases, however, it may be difficult to estimate angular velocity only from images taken at one time (e.g. due to poor illumination of the object 50, or the object 50 being too distant from the cameras). The inventors have realised that, by determining the angular velocity at a plurality of times (from three or more images captured at different times), the reliability and / or accuracy of an estimate made for a given time can be improved, as follows. The angular velocity of the ball can be affected by air resistance and by bouncing. However, in general, the angular velocity of the ball is unlikely to change significantly in a short period of time. To import this observation into the present angular velocity estimation method, the method 100 uses a modified cost function, which takes account of a second angular velocity wt2 of the object 50 for a time t2 different to t0. Where the second angular velocity wt2is a known value, the modified cost function P’ is , is referred to as a smoothing term. For example, the angular velocity wta for a time ta may be estimated using the non-modified cost function with images captured at times taand tb. Then, to estimate the angular velocity wtbfor a time tbfrom the image captured at time tband a third image captured at a later time tc, the previously estimated angular velocity wta may be treated as a known value wt2 and used in the modified cost function shown above, with the photogrammetric term using the pair of images captured at times tband tc. Since the modified cost function includes a term that increases with |wta - wtb|, the determined value of wtb is more likely to be closer to than if the smoothing term is absent. (Where the angular velocities wtaand wtbare represented as matrices R(θ) (representing a rotation transformation, around the candidate rotation axis of the object 50, by an angle θ, equal to e.g. the angular velocity of the object 50 multiplied by 1 second), the term |wta- wtb| can be calculated as acos((tr(R(θta)T R(θtb))-1) / 2).) The function is a positive, monotonically decreasing function of (t2– t0). This function takes account of the lack of knowledge of the torque acting on the object 50 between t0 and t2: the further apart t0 and t2 are, the more likely it is that the true angular velocity of the object 50 has in fact changed significantly, and hence the less the “smoothing” should be taken account of. The function may be, for example, an exponential decay function i.e. exp(-|t2 – t0|), a Gaussian function exp(-(t2 – t0)2), or a function defined as 1 / (1 +|^^− ^^|). The above form of the modified cost function assumes that the second angular velocity wt2of the object 50 is known. However, even when no angular velocity of the object 50 at any time is known, the cost function can be similarly modified, if angular velocities are to be calculated for two or more times. The generalised form of the cost function is where wt0,…,wtn are the angular velocities of the object 50 at respective times t0, t1,…,tn, I(x0{i0}) is the intensity of the one or more pixels in region x0in a first image i0 (e.g. first image 40) of an image pair i, and I(x1{i1}) is the intensity of the one or more pixels in corresponding region x1in a second image i1(e.g. second image 41) of the image pair i. Note that the images in image pairs i could be taken at any of the times t0, (or even at times between these times), and some of these times might be covered by multiple (e.g.2 or more) images. By calculating two angular velocities for respective times using a joint cost function with a smoothing term, as described above, a stable and accurate estimate of the angular velocity can be obtained for each of the times. The above embodiments are to be understood as illustrative examples of the invention. It is to be understood that any feature described in relation to any one embodiment may be used alone, or in combination with other features described, and may also be used in combination with one or more features of any other of the embodiments, or any combination of any other of the embodiments. Furthermore, equivalents and modifications not described above may also be employed without departing from the scope of the invention, which is defined in the accompanying claims.

Claims

CLAIMS 1. A method for estimating an orientation change characteristic of an object from a plurality of images depicting the object captured at respective different times, the method comprising: for each of a plurality of candidate values of the orientation change characteristic for the object: identifying a first region in a first image of the plurality of images and a second region in a second image of the plurality of images, wherein the first region and the second region would represent a same portion of a surface of the object if the orientation change characteristic was equal to the candidate value; and determining, based on pixel intensity values of one or more pixels in the first region, pixel intensity values of one or more pixels in the second region, and an illumination parameter set comprising a respective value for each of one or more illumination parameters, the illumination parameter set representing illumination conditions of the object in the plurality of images, a value of a cost function; and estimating the orientation change characteristic based on the values of the cost function.

2. The method of claim 1, wherein determining the value of the cost function comprises: determining, based on the illumination parameter set, at least one illumination- adjusted difference between the pixel intensity values of the one or more pixels in the first region and the pixel intensity values of the one or more pixels in the second region; and determining the value of the cost function based on the illumination-adjusted difference.

3. The method of claim 1 or claim 2, wherein the illumination parameter set is a first candidate illumination parameter set of a plurality of candidate illumination parameter sets, and the method comprises:for each of the plurality of candidate values of the orientation change characteristic for the object: for each of the plurality of candidate illumination parameter sets, determining, based on the pixel intensity values of the one or more pixels in the first region, the pixel intensity values of the one or more pixels in the second region, and the candidate illumination parameter set, an illumination-specific value of the cost function; and estimating the orientation change characteristic based on the illumination- specific values of the cost function.

4. The method of claim 3, comprising determining the illumination parameter set representing the illumination conditions of the object in the plurality of images by selecting a candidate illumination parameter set, from among the plurality of candidate illumination parameter sets, based on the values of the cost function.

5. The method of any one of claim 1 to claim 4, wherein estimating the orientation change characteristic comprises selecting a candidate value of the orientation change characteristic, from among the plurality of candidate values, that has a lowest value of the cost function among the values of the cost function.

6. The method of any one of claim 1 to claim 5, wherein the illumination parameters comprise at least one of: a location of a light source relative to the object, a direction of illumination of the object, an intensity of the illumination, a colour of the illumination, a specular reflection property of the object, a diffuse reflection property of the object, or an ambient reflection property of the object.

7. The method of any one of claim 1 to claim 6, wherein the orientation change characteristic is a first angular velocity for a first time, and estimating the first angular velocity comprises, for each of the plurality of candidate values of the first angular velocity: determining the value of the cost function further based on a difference between the candidate value of the first angular velocity and a second angular velocity of theobject for a second, different, time, so as to reduce a difference between the estimated first angular velocity and the second angular velocity.

8. The method of claim 7, wherein determining the value of the cost function comprises weighting the difference between the candidate value of the first angular velocity and the second angular velocity according to a difference between the first time and the second time.

9. The method of any one of claim 1 to claim 8, wherein identifying the second region comprises: determining a first position of the portion of the surface of the object in 3D space, from the first image; determining, based on the first position and the candidate value of the orientation change characteristic, a second position of the portion of the surface of the object in the 3D space, the second position representing a position of the portion of the surface of the object in the 3D space at a time of capture of the second image; and determining, based on the second position of the portion of the surface of the object in the 3D space, the second region in the second image.

10. The method of claim 9, wherein determining the first position comprises determining the first position based on a known size of the object and based on a measured size of the object in the first image, and / or identifying the second position comprises identifying the second position based on the known size of the object and based on a measured size of the object in the second image.

11. The method of any one of claim 1 to claim 10, wherein the orientation change characteristic comprises an angular speed, at least one component of angular velocity, or a direction of angular velocity, of the object.

12. The method of any one of claim 1 to claim 11, wherein the object is a sports ball.

13. The method of any one of claim 1 to claim 12, comprising: determining, based on a time difference between the first image and the second image, and on a maximum expected angular speed of the object, that a maximum expected angular displacement of the object between the first image and the second image is lower than a threshold, wherein identifying the first region in the first image and the second region in the second image is performed in response to determining that the maximum expected angular displacement is lower than the threshold.

14. A processing system configured to perform the method of any one of claim 1 to claim 13.

15. A system comprising: the processing system of claim 14, and one or more cameras, the one or more cameras being configured to provide the plurality of images to the processing system.

16. A non-transitory computer readable medium comprising instructions which, when executed by a processor, cause the processor to perform the method of any one of claim 1 to claim 13.

Images