Daniel Krpelik

Milan Jaros

Marta Jarosova

Petr Strakos

Tereza Buresova

Alena Vasatova

Tomas Karasek
email: tomas.karasek@vsb.cz

Detection of Orbital Floor Fractures by Principal Component Analysis

Principal component analysis (PCA) is a statistical method based on orthogonal transformation, which is used to convert possibly correlated datasets into linearly uncorrelated variables called principal components. PCA is one of the simplest methods based on the eigenvector analysis. This method is widely used in many fields, such as signal processing, quality control or mechanical engineering. In this paper, we present the use of PCA in area of medical image processing. In the medical image processing with subsequent reconstruction of 3D models, data from sources such as Computed Tomography (CT) or Magnetic Resonance Imagining (MRI) are used. Series of images representing axial slices of human body are stored in Digital Imaging and Communications in Medicine (DICOM) format. Physical properties of different body tissues are characterized by different shades of grey of each pixel correlated to the tissue density. Properties of each pixel are then used in image segmentation and subsequent creation of 3D model of human organs. Image segmentation splits digital image into regions with similar properties which are later used to create 3D model. In many cases accurate detections of edges of such objects are necessary. This could be for example the case of a tumour or orbital fracture identification. In this paper, identification of the orbital fracture using PCA method is presented as an example of application of the method in the area of medical image processing.

Introduction

Despite CT and MRI are widely used by medical doctors for diagnostic purposes their potential is still not fully utilized. Those modalities provide 2D images of 3D objects and although medical doctors learn how to read and interpret them, they provide only partial information. Extraction of additional information such as size, area or volume of desired object has to be done manually and it is labour-intensive. Information about volume is crucial for example for surgeons, who are planning liver resections where remaining volume of liver is matter of life or death. Volume of upper respiratory tract serves in diagnoses of sleep apnea. Size of orbital floor fracture is important when doctors are making decision whether patient should undergo surgery or conservative treatment. Those are just examples of additional information requested by doctors which cannot be obtained automatically from CT or MRI images.

Image processing could be an option how to reduce labour intensity, shorten the time of diagnostic processes and increase fidelity of the results. Above mentioned quantities are measured on 3D models which are created usually in several following steps. First, image noise is reduced by image filtering. Than image segmentation consisting of localization of objects and its boundaries is performed. Last step is 3D surface reconstruction.

If the 3D model is created by the mentioned steps, further detection of significant features like edges can be highly important. These features can localize major changes in overall shape of the model. While working with medical data, localization of such changes can lead, e.g., to direct detection of fractures in bones or tumours in soft tissues.

PCA method is used here to detect the points on the edge of the fracture. The edges are identified by calculating eigenvalues and eigenvectors of covariance matrix and establishing ratios between the calculated eigenvalues. In this paper PCA method is applied for detection of orbital floor fracture.

This paper is organized as follows. In Section 2 we describe example of identification of orbital floor fracture. In Section 3 PCA method is described followed by chapter 4 where object identification is explained. In chapter 5 details about fracture area reconstruction can be found. Chapter 6 concludes whole paper.

Model example

In this paper PCA is used for identification of orbital floor fracture from 3D virtual model of the human orbit. There is a very thin bone on the bottom side of the human eye orbit. Since this bone is the most fragile part of the orbit it is predisposed to fractures (see Fig. 1, Fig. 2). Fracture of this bone can be determined from CT scans, but to determine the extension of the fracture, only rough approximate methods have been used so far. These methods rely on manual determination of most distant points on fracture boundary. The fracture size can then be assessed by calculating the ratio of distance between two points of fracture boundary, to the distance between another two points, representing orbital floor boundary. All point are selected in one CT image where the fracture diameter is maximal [START_REF] Ploder | Evaluation of computer-based area and volume measurement from coronal computed tomography scans in isolated blowout fractures of the orbital floor[END_REF]. Another way is to select points denoting fracture boundary in two perpendicular planes and approximate the fracture by geometrical primitives (e.g. ellipse). Area of such primitives can be determined from known equations [START_REF] Goggin | Simple Computed Tomography-Based Calculations of Orbital Floor Fracture Defect Size Are Not Sufficiently Accurate for Clinical Use[END_REF]. In this approach CT image with maximal fracture diameter is also used. Although both above mentioned methods were shown to be inaccurate, they are widely used as they require little human expert intervention and less time compared to purely manual segmentation. Manual segmentation requires manual selection of all points denoting the fracture on each individual CT image followed by area calculation from selected voxels. Our implemented method is based on semi-automatic segmentation and keeps the human expert intervention minimal.

Covariance matrix and principal component analysis

First, let us introduce the basic terms which will be used later in this paper. To distinguish, whether particular object is 1D, 2D or 3D we need to know the relationship between points representing this object. Let P ⊂ R 3 be the set of all data points. For each point p ∈ P set of neighbouring nodes could be found as

N p = q | p, q ∈ P, |pq| = (p x -q x ) 2 + (p y -q y ) 2 + (p z -q z ) 2 < d 3D . (1)
The example of set of points for which covariance matrix will be created is depicted on Fig. 3.

a) p Np l v 1 1 l v 2 2 l v 3 3 b) p Np l v 1 1 l v 2 2 l v 3 3 N p c) p l v 1 1 l v 2 2 l v 3 3 Fig. 3. Point distribution: a) linear, b) planar, c) spherical
Recognition of geometric patterns in an image by the invariant moment was first introduced by 1962 Hu [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF]. In this work two-dimensional moment was used. Later, the invariant moment was extended to the third dimension by [START_REF] Maas | Two algorithms for extracting building models from raw laser altimetry data[END_REF] and [START_REF] Gross | Extraction of lines from laser point clouds[END_REF].

Let us consider the set N p to be solid in such a way that B ⊂ R 3 . Then (i + j + k)-th moment is defined by the equation

m ijk = B x i y j z k f (x, y, z) dxdydz, (2)
where i, j, k ∈ Z + 0 and f (x, y, z) is the mass density function of solid B. Let f (x, y, z) = 1 for ∀p ∈ P .

If we set (i + j + k) ≤ 2, coordinates of the centre of gravity p of the solid B could be defined as

p = x = m 100 m 000 , ȳ = m 010 m 000 , z = m 001 m 000 , (3)
where m 000 is the mass of the solid B and m 100 , m 010 , m 001 are the static moments with respect to the coordinate planes yz, xz, xy. The (i + j + k)-th centralized moment of the solid B is defined by

µ ijk = B (x -x) i (y -ȳ) j (z -z) k f (x, y, z) dxdydz. (4)
To obtain central moment invariant to its size we have to normalize it

µ ijk = µ ijk µ 000 = B (x -x) i (y -ȳ) j (z -z) k f (x, y, z) dxdydz B (x -x) 0 (y -ȳ) 0 (z -z) 0 f (x, y, z) dxdydz . ( 5
)
After the numerical approximation of the normalized centralized moment we get

µ ijk = |Np| n=1 (x n -x) i (y n -ȳ) j (z n -z) k |Np| n=1 (x n -x) 0 (y n -ȳ) 0 (z n -z) 0 , ( 6
)
µ ijk = 1 |N p | |Np| n=1 (x n -x) i (y n -ȳ) j (z n -z) k . (7)
Finally, covariance matrix for each point p ∈ P and its neighbours could be assembled as

C p =   µ 200 µ 110 µ 101 µ 110 µ 020 µ 011 µ 101 µ 011 µ 002   . (8)
The matrix C p could be simplified by substituting the moments

C p = 1 |N p | |Np| n=1   (x n -x) 2 (x n -x)(y n -ȳ) (x n -x)(z n -z) (x n -x)(y n -ȳ) (y n -ȳ) 2 (y n -ȳ)(z n -z) (x n -x)(z n -z) (y n -ȳ)(z n -z) (z n -z) 2   , (9)
C p = 1 |N p | |Np| n=1     x n - x y n - ȳ z n - z   x n -x y n -ȳ z n -z   , (10)
C p = 1 |N p | |Np| n=1 (q n -p) T (q n -p), (11)
where q ∈ N p and p is the center of gravity of the set N p . For purpose of implementation, equation (11) could be rewritten to the form of matrix multiplication

Q p =        x 1 y 1 z 1 x 2 y 2 z 2 x 3 y 3 z 3 . . . . . . . . . x |Np| y |Np| z |Np|        , J =        1 1 . . . 1 1 1 . . . 1 1 1 . . . 1 . . . . . . . . . 1 1 . . . 1        , (12)
S p = Q p -J T Q p 1 |N p | , (13)
C p = S T p S p , (14)
where Q p is the matrix composed of coordinates of a set N p , matrix J is square matrix of all ones and the matrix multiplication J T Q p 1 |Np| represents the center of the gravity of the set N p .

The covariance matrix C p contains only real numbers, it is symmetric positive semi-definite and its dimension is n = 3. It means that all its eigenvalues are real with a complete set of the orthonormal eigenvectors. If this is a case we could say that there exists an orthogonal matrix Q and a diagonal matrix D such that

C p = QDQ T . ( 15
)
The diagonal elements of the matrix D are the eigenvalues λ 1 ≥ λ 2 ≥ λ 3 and columns of the matrix Q are the orthonormal eigenvectors v 1 , v 2 , v 3 of the matrix C p

C p = v T 1 v T 2 v T 3   λ 1 0 0 0 λ 2 0 0 0 λ 3     v 1 v 2 v 3   . ( 16
)

Object recognition

An auxiliary 3D model is acquired by prior segmentation of CT images. Surface of this model is then converted into cloud of points. Following task is to select points which lie on the fracture boundary. This is solved by the classification of neighbourhood of each single point in the cloud. Auxiliary model surface is a 2D object in 3D space. Points of the fracture boundary are determined based on violation of planarity in their neighbourhoods. Object recognition, respectively classification whether object is 1D, 2D or 3D is based on properties of eigenvalues λ 1 ≥ λ 2 ≥ λ 3 . Eigenvalues are obtained by spectral decomposition of matrix C p in equation ( 14) and they will represent the shape of the set N p by the ellipsoid [START_REF] Gross | Extraction of lines from laser point clouds[END_REF][START_REF] Demantke | Dimensionality based scale selection in 3D lidar point clouds[END_REF]. Ellipsoid representing linear, planar or spheric behaviour of neighbours is depicted in Fig. 3. If λ 1 >> λ 2 , λ 3 ∼ = 0, the point p is marked as linear. If λ 1 , λ 2 >> λ 3 ∼ = 0, the point p is assigned to the planar group. If λ 1 ∼ = λ 2 ∼ = λ 3 , the point p belongs to spheric group.

For the detection of a fracture we divide the points into two groups, P planar and P boundary using the criteria

a planar = λ 3 d λ d , (17)
P planar = {p|p ∈ P : a planar ≤ } , (18)
P boundary = {p|p ∈ P : a planar > } , (19)
where d λ d is the normalization coefficient and is the tolerance parameter. If a planar > , the point is on the edge and it is assigned to the group P boundary .

Otherwise the point is on the plane and it is assigned to the group P planar .

Reconstruction of fracture area

Fracture is in fact the absence of bone and as such it cannot be localized based solely on pixel/voxel intensity values but rather on its relative position to the other parts of human body. The course of work involves prior segmentation and reconstruction of either undamaged bone area or eye orbit interior as sought fracture boundary lies on the intersection of their surfaces. The eye orbit interior is composed of soft tissue and in many cases leaks through the fracture into air cavity bellow, creating a droplet-like structure. The fracture boundary itself is then searched as a distinctive subsection of selected auxiliary object surface, mainly determined by relatively extreme curvature in its neighbourhood.

This curvature is then exploited in two ways. Firstly, in case of both manual and automated point selection, this curvature is used for reconstruction of fracture boundary. Working hypothesis is that the shortest Jordan curve on the object surface passing through selected points is a good approximation of fracture boundary. Secondly, since the curvature is a distinctive feature, it may be used to extract boundary points via automated segmentation of auxiliary model vertices.

In applied method, we are reconstructing the fracture boundary. Fracture interior itself is impossible to determine, because we cannot evaluate, which specific surface shape is the correct one. Fracture interior is reconstructed heuristically as a piecewise planar surface bounded by found fracture boundary. Alternative methods could prefer smooth surfaces and reconstruct fracture interior for example via solving Laplace problem with Dirichlet boundary.

Proposed method for Jordan curve localization assume that selected points are vertices of auxiliary object mesh and are distributed sparsely around sought boundary. The points in the case of manual point selection are selected on the CT scan images and might not correspond with any vertex in the auxiliary object mesh. This is solved simply by selecting the closest vertex in the mesh for each manually selected point. In case of automated point selection we usually acquire a dense cloud of points, which denies requirement on sparsity. This might be solved automatically by, e.g., the k-means clustering, in which case we select output k-means centroids positions and treat those as a compression of initially dense cloud of points.

The following boundary reconstruction consists of two steps. First we construct a rough approximation of boundary consisting only of selected points and edges, and directly connecting them across Euclidean space. In terms of creating boundary, these points have to be properly ordered. This is a special case of travelling salesman problem for Euclidean metric space, which is generally NP-hard. Although a good polynomial approximation is known, we, for the sake of simplicity, use a greedy algorithm.

Our greedy algorithm starts with arbitrary triangle (a cycle with 3 vertices) and iteratively lengthens the cycle by inserting new point in between a pair of consecutive points from previous iteration. The pair is selected in such a way that the resulting new cycle is the shortest one among those considered (length being measured as a sum of distances between all consecutive points in the cycle). This algorithm is polynomial and given our assumptions, it yields good approximation of sought shortest cycle.

As a next step we refine found boundary by replacing edges in between consecutive points by the shortest paths between them on the auxiliary object mesh. This can be done effectively via Dijkstra algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] or some of its approximating variants. To ensure properties of Jordan curve a simple postprocessing consisting of removing the inner cycles, which may occur because of independent runs of Dijkstra algorithms, is carried out.

Resulting reconstructed fracture area is visualized on Fig. 5 and Fig. 6.

Conclusion

A method for edge recognition applied to detect orbital floor fracture was examined in this paper. The method is based on utilization of auxiliary 3D models obtained by segmentation of CT scan images. The problem of fracture reconstruction was converted to a problem of segmentation of auxiliary model mesh vertices and edges.

Proposed method requires selection of specific vertices on auxiliary model mesh. This selection was carried out both manually and automatically, via discriminating features based on PCA, and the results have been compared visually. Once the points were selected, algorithm chose appropriate subset of edges from auxiliary model mesh to reconstruct fracture boundary. These edges were selected so that they correspond with the shortest paths between selected points. This was carried out by Dijkstra algorithm.

Although the throughout validation of the size of orbital floor fracture based on semi-automatic segmentation is still pending, preliminary results based on expert based validation of fracture areas projected backward into CT images show promising results. PCA based point selection provides satisfactory results and thus it deserves further development. For example, the used discriminative feature, the relative magnitude of the smallest eigenvalue, is unable to fully capture all the parts of the fracture boundary. Finding better discriminative feature which could overcome this problem will be our future work.

Fig. 1 .

 1 Fig. 1. Model example: Left image depicts eye area denoted by yellow rectangle on CT scan image of a skull. Right image shows boundary (pink colour) of inner eye region identified by segmentation algorithm.

Fig. 2 .

 2 Fig. 2. 3D reconstruction of orbital floor. Hole on the bottom side represents sought fracture.

Fig. 4 .

 4 Fig. 4. Conversion of voxels to point cloud (left) selection of cloud subset highlighted in green color (middle) and result of automatic boundary points extraction (right)

Fig. 5 .

 5 Fig. 5. Fracture extraction -backward projection: manual point selection (pink), automatic point selection (yellow)

Acknowledgement

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project "IT4Innovations excellence in science -LQ1602" and from the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National Supercomputing Center -LM2015070". All CT images used in this paper were provided by The University Hospital in Ostrava.


Lukáš Vojáček
email: lukas.vojacek@vsb.cz

Pavla Dráždilová
email: pavla.drazdilova@vsb.cz

Jiří Dvorský
email: jiri.dvorsky@vsb.cz

Optimization of combining of Self Organizing Maps and Growing Neural Gas

Keywords: self organizing maps, growing neural gas, high-dimensional dataset, high performance computing

The paper deals with the issue of high dimensional data clustering. One possible way to cluster this kind of data is based on Artificial Neural Networks (ANN) such as Growing Neural Gas (GNG) or Self Organizing Maps (SOM). Parallel modification, Growing Neural Gas with pre-processing by Self Organizing Maps, and its implementation on the HPC cluster is presented in the paper. Some experimental results are also presented. We focus on effective preprocessing for GNG. The clustering is realized on the output layer of SOM and the data for GNG are distributed into parallel processes.

Introduction

Recently, several topics and issues have emerged such as availability of big and/or high-dimensional data, petascale HPC systems aiming towards exascale supercomputers and need of processing this kind of data. Large high-dimensional data collections are commonly available in areas like medicine, biology, information retrieval, web analyze, social network analyze, image processing, financial transaction analysis and many others. To process such kind of data unsupervised learning algorithms, such as Self Organizing Maps (SOM) or Growing Neural Gas (GNG), are usually used. Various aspects of parallel implementation of these algorithms on HPC were studied e.g. [START_REF] Vojáček | Parallel hybrid som learning on high dimensional sparse data[END_REF][START_REF] Vojáček | Combination of self organizing maps and growing neural gas[END_REF].

The one of still preserving issues is efficient utilization of the computation resources. When speaking on SOM or GNG learning algorithms there are two challenges. The first one is fast computation of similarity in high-dimensional space and the second one is ideally uniform distribution of computation load among individual CPU cores.

Parallel implementation usually allocates one CPU core to group of neurons, evaluate similarity of these neurons with given input vector, find local best matching neuron and then using some form of communication to find a global best matching neuron in the whole neural network. The CPU cores are allocated regularly, using some pattern [START_REF] Vojáček | Combination of self organizing maps and growing neural gas[END_REF] regardless of input vectors distribution on neurons i.e. CPU cores causing a bottleneck in the parallel learning algorithm. To reduce the bottleneck input vectors preprocessing is done using small SOM and clustering algorithm. This allows us to improve distribution of neurons over CPU cores and subsequently speed-up the learning algorithm itself.

The paper is organized as follows. The Section 2 briefly describes used neural networks: SOM and GNG. Input vectors preprocessing algorithm is provided in Sect. 3. Experimental results are given in Sect. 4.

Artificial Neural Networks

In this section we will describe two types of neural networks, the first is Self Organizing Maps and the second is Growing Neural Gas and then we present a combination of SOM and GNG.

Self Organizing Maps

Self Organizing Maps (SOMs), also known as Kohonen maps, were proposed by Teuvo Kohonen in 1982 [START_REF] Kohonen | Self-Organization and Associative Memory[END_REF]. SOM is a kind of artificial neural network that is trained by unsupervised learning. Using SOM, the input space of training samples can be represented in a lower-dimensional (often two-dimensional) space [START_REF] Kohonen | Self Organizing Maps[END_REF], called a map. Such a model is efficient in structure visualization due to its feature of topological preservation using a neighbourhood function.

SOM consists of two layers of neurons: an input layer that receives and transmits the input information, and an output layer, the map that represents the output characteristics. The output layer is commonly organized as a twodimensional rectangular grid of nodes, where each node corresponds to one neuron. Both layers are feed-forward connected; each neuron in the input layer is connected to each neuron in the output layer. A real number, or weight, is assigned to each of these connections.

Growing Neural Gas

The representation of Growing Neural Gas is an undirected graph which need not be connected. Generally, there are no restrictions to the topology. The graph is generated and continuously updated by competitive Hebbian Learning [START_REF] Martinetz | Competitive hebbian learning rule forms perfectly topology preserving maps[END_REF][START_REF] Prudent | An incremental growing neural gas learns topologies[END_REF]. According to the pre-set conditions, new neurons are automatically added and connections between neurons are subject to time and can be removed. GNG can be used for vector quantization by finding the code-vectors in clusters [START_REF] Holmström | Growing Neural Gas Experiments with GNG, GNG with Utility and Supervised GNG[END_REF], image compression, disease diagnosis.

GNG works by modifying the graph, where the operations are the addition and removal of neurons and edges between neurons.

To understand the functioning of GNG, it is necessary to define the learning algorithm. The algorithm published in [START_REF] Vojáček | Combination of self organizing maps and growing neural gas[END_REF] is based on the original algorithm [START_REF] Fritzke | A growing neural gas network learns topologies[END_REF][START_REF] Holmström | Growing Neural Gas Experiments with GNG, GNG with Utility and Supervised GNG[END_REF], but it is modified for better continuity in the SOM algorithm.

Remark The notation used in the paper is briefly listed in Table 1.

Combination of SOM and GNG

In our previous paper [START_REF] Vojáček | Combination of self organizing maps and growing neural gas[END_REF], we focused on a combination of SOM and GNG, where the basic idea was to pre-process the input data by SOM, as a result of which there are clusters of similar data. Subsequently, we created the same number of GNG network as clusters, and assigned each cluster to one GNG. Each GNG creates its own neural map and after the learning process is finished, the results are merged. The entire description above can be summarized as follows: Help speeding up computation parallelization is shown in Fig. 1 where the top layer of parallelization (SOM) is described in a previous paper [START_REF] Vojáček | Parallel hybrid som learning on high dimensional sparse data[END_REF]. In this chapter we will describe an improved method focusing on the creating clusters of input data and optimization of used resources.

To improve the efficiency of parallelization is needed to clusters of input data for GNG network contained approximately the same number of input data. Based on the fact that when the GNG assigned more input vectors, the calculation takes longer. In the past to create clusters we used to spanning tree algorithm [START_REF] Vojáček | Combination of self organizing maps and growing neural gas[END_REF] which, however, does not reflect the number of input vectors in clusters. To obtain clusters of neurons we now use two different clustering algorithms. The first algorithm is agglomerative hierarchical clustering methods We proposed an optimization for the calculation of GNG networks, which aims to optimize the maximum utilization of the allocated resources, but on condition that the computing time must be similar. The principle of optimization is based on the idea that individual computing resources will count more GNG networks than only one -as it has until now. In Algorithm 2, the overall functionality is described. Here it is necessary to mention two facts regarding point 4. Firstly, it is a variation of a known problem Subset sum [START_REF] Li | On the subset sum problem over finite fields[END_REF], which is an NP complete problem [START_REF] Carlson | The Millennium Prize Problems[END_REF]. But at the beginning it was defined that the size of the output map of SOM is small and therefore the maximum possible number of clusters is also small, and therefore negligible. And secondly the maximum limit is reduced by 10% on the grounds that to work with each GNG network time for I/O operations must be added. where |Ni| denotes number of input vectors mapped to neuron Ni. 2. Find a K ∈ {5, 6, 7, 8} for which value N k F -1 k 2 M is minimal. So that SOM (K × K) has the smallest dispersion of |Ni|. Furthermore, we will work with SOM (K × K). For sake of simplicity the index K will be omitted from N K F in following text.

V16 V22 V17 V21 V6 V12 V7 V11 V1 V2 V13 V8 V14 V3 V9 V10 V4 V5 V15 V20 V25 V18 V24 V19 V23 0 2 4 6 8 V10 V4 V5 V15 V3 V9 V25 V20 V24 V8 V14 V13 V19 V1 V2 V6 V12 V7 V11 V23 V18 V22 V16 V17

Experiments

We will describe different datasets and we will provide experiments with datasets in this section.

Experimental Datasets and Hardware

Two datasets were used in the experiments. The first dataset was commonly used in Information Retrieval -Medlars. The second one was the test data for the elementary benchmark for clustering algorithms [START_REF] Ultsch | Clustering with SOM: U*C[END_REF].

Weblogs Dataset To test the learning algorithm effectiveness on high dimensional datasets, a Weblogs dataset was used. The Weblogs dataset contained web logs from an Apache server. The dataset contained records of two month's requested activities (HTTP requests) to the NASA Kennedy Space Center WWW server in Florida3 . The standard data preprocessing methods were applied to the obtained dataset. The records from search engines and spiders were removed, and only the web site browsing was left (without download of pictures and icons, stylesheets, scripts etc.). The final dataset (input vector space) was of a dimension 90,060 and consisted of 54,961 input vectors. For a detailed description, see our previous work [START_REF] Slaninová | Web site community analysis based on suffix tree and clustering algorithm[END_REF], where a web sites community behaviour was analyzed.

Medlars Dataset

The Medlars dataset consisted of 1,033 English abstracts from medical science4 . The 8,567 distinct terms were extracted from the Medlars dataset. Each term represents a potential dimension in the input vector space. The term's level of significance (weight) in a particular document represents a value of the component of the input vector. Finally, the input vector space has a dimension of 8,707, and 1,033 input vectors were extracted from the dataset.

Experimental Hardware

The experiments were performed on a Linux HPC cluster, named Anselm, with 209 computing nodes, where each node had 16 processors with 64 GB of memory. Processors in the nodes were Intel Sandy Bridge E5-2665. Compute network is InfiniBand QDR, fully non-blocking, fattree. Detailed information about hardware is possible to find on the web site of Anselm HPC cluster 5 .

First Part of the Experiment

The first part of the experiments was oriented towards the comparison of quality of clustering by the agglomerative hierarchical clustering and PAM Clustering. The used dataset was Weblogs. All the experiments were carried out for 20 epochs; the random initial values of neuron weights in the first epoch were always set to the same values. The tests were performed for SOM with rectangular shape 5 × 5 neurons, 6 × 6 neurons, 7 × 7 neurons and 8 × 8 neurons. The metrics used for determining the quality of clustering is Average Silhouette Index (ASI). The achieved quality of clustering for agglomerative hierarchical clustering is presented in Fig. 3 and for PAM Clustering is presented in Fig. 4.

Second Part of the Experiment

The second part of the experiments was oriented towards comparing the time efficiency of algorithms PAM and AVL. The GNG parameters are as follows γ = 100, e w = 0.05, e n = 0.006, α = 0.5, β = 0.0005, a max = 160, N max = 1221, T = 200. The used dataset was Weblogs. Dimensions of SOM are 5 × 5, 6 × 6, 7 × 7 and 8 × 8. Number of cores are 32 for each group and the groups are computed sequentially. The tests were performed for SOM with a rectangular shape of 5 × 5 neurons. The achieved speedup is presented in Table 4 and Table 5. Both samples have similar speedup but fundamentally differ in efficiency; in the fastest case without optimization the efficiency is 0.02, but with optimization it is 0.23.

Conclusion

In this paper, we presented experiments with preprocessing of data in output layer of SOM. Different clustering algorithm was used and different quality of clusters were obtained. The global input data preprocessed on this SOM was used as an input for GNG neural network. This approach allows us almost uniformly distribute data on computation cores efficiently utilize them. The achieved speedup is also very good. In future work we intend to focus on sparse date, and improved acceleration and use Xeon phi for better speed-up.

Fig. 1 .

 1 Fig. 1. Parallel algorithm

Fig. 2 .Algorithm 1 :

 21 Fig. 2. Dendrograms of hierarchical clustering (Average-linkage and Ward's method) of output layer of SOM 5x5

3 .

 3 Divide the cluster C i with the fattest neuron NF to two clusters C N F and Ci where cluster C N F contain only one neuron NF and cluster Ci = C i \{NF }. 4. For others clusters C j verify: if number of input data |C j | in cluster C j |C j | > |NF |, then split C j into two clusters Cj,1 and Cj,2 where cluster Cj,1 contain fattest neuron from C j and cluster Cj,2 = C j \Cj,1. 5. Return new set of cluster.

Algorithm 2 :

 2 Algorithm to optimize the utilization of computing resources 1. From the set of clusters C = {C1, . . . , Cn}; take cluster C N F with the greatest number of input vectors. 2. Remove C N F from C. 3. Assign cluster C N F to unused computing resource. 4. Select the set of clusters from set C wherein the sum of their input vectors is approaching number of input vectors in 0.9|C N F | 5. Remove used clusters from the set C. 6. Assign selected clusters to the unused computing resource. 7. If C is not empty return to step 4.

Fig. 4 .

 4 Fig. 4. Quality of PAM clustering of output layer of SOM k×k

Table 1 .

 1 Notation used in the paper

	Symbol Description
	M	Number of input vectors
	n	Dimension of input vectors, number of input neurons, dimension of weight
		vectors in GNG output layer neurons
	N	Current number of neurons in GNG output layer
	Nmax Maximum allowed number of neurons in GNG output layer
	ni	i-th input neuron, i = 1, 2, . . . , n
	Ni	i-th output neuron, i = 1, 2, . . . , N
	NF	the fattest neuron
	T	Number of epochs
	lc 1	Learning factor of BM U1
	lnc

1 Learning factor of BM U1 neighbours ei Local error of output neuron Ni, i = 1, 2, . . . , N Ci ith cluster witch contains similar input vectors mi Number of input vectors in cluster Ci Zi Centroid of cluster Ci α Error ei reduction factor β Neuron error reduction factor γ Interval of input patterns to add a new neuron amax Maximum edge age

Table 2 .

 2 Computing time with respect to number of cores, standard GNG algorithm,

	dataset Medlars
	SOM		Numbers of vectors		Time [mm:ss]
	Dimension Time [s] Max Second max 4* PAM 10 PAM 9* AVL 4 AVL
	5 × 5	13	6852	1636	14:19	14:37	14:44	14:22
					6* PAM 10 PAM 11* AVL 6 AVL
	6 × 6	18	5918	1900	13:50	14:12	14:11	14:02
					10* PAM 6 PAM 20* AVL 6 AVL
	7 × 7	27	5607	1795	13:57	13:47	14:26	13:48
					9* PAM	6 PAM 20* AVL 6 AVL
	8 × 8	35	5192	1437	14:03	13:51	14:41	13:37
	5 https://support.it4i.cz/docs/anselm-cluster-documentation/hardware-overview

Table 3 .

 3 The quality of algorithms (the lower is better)

	4* PAM 10 PAM 9* AVL 4 AVL
	1.24	1.23	1.24	1.2
	6* PAM 10 PAM 11* AVL 6 AVL
	1.22	1.24	1.22	1.26
	10* PAM 6 PAM 20* AVL 6 AVL
	1.14	1.38	1.4	1.18
	9* PAM	6 PAM 20* AVL 6 AVL
	1.24	1.04	1.3	1.18

Table 4 .

 4 Combination

				SOM and
	GNG's
	Cores	Time	Speedup Efficiency
	1 02:28:46	1.00	1.00
	32 00:18:01	8.26	0.26
	64 00:17:51	8.33	0.13
	128 00:16:37	8.95	0.07
	256 00:10:42	13.90	0.05
	512 00:07:27	19.97	0.04
	1024 00:06:37	22.48	0.02

Table 5 .

 5 Combination SOM and GNG's and resource optimization

	Cores	Time	Speedup Efficiency
	1 02:28:46	1.00	1.00
	32 00:18:02	8.26	0.26
	64 00:17:48	8.33	0.13
	96 00:06:39	22.37	0.23
	-	-	-	-
	-	-	-	-
	-	-	-	-

The collection can be downloaded from http://ita.ee.lbl.gov/html/contrib/NASA-HTTP.html.

[START_REF] Kohonen | Self-Organization and Associative Memory[END_REF] The collection can be downloaded from ftp://ftp.cs.cornell.edu/pub/smart. The total size of the dataset is approximately 1.03 MB.

Acknowledgment

This work was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National Supercomputing Center LM2015070 " and partially supported by Grant of SGS No. SP2016/68, V ŠB -Technical University of Ostrava, Czech Republic.


Piotr Jurewicz

Wojciech Rafaj Lowicz

Jacek Reiner
email: jacek.reiner@pwr.edu.pl

Ewaryst Rafaj Lowicz
email: ewaryst.rafajlowicz@pwr.edu.pl

Simulations for Tuning a Laser Power Control System of the Cladding Process

Keywords: control system, simulations, IR camera, laser cladding

Our aim is present the methodology of simulations for repetitive processes and tuning control systems for them in the presence of noise. This methodology is applied for tuning a laser power control system of the cladding process. Even the simplest model of this process is nonlinear, making analytical tuning rather difficult. The proposed approach allows us to select quickly the structure of the control system and to optimize its parameters. Preliminary comparisons with experimental results on a robot-based laser cladding systems are also reported. These comparisons are based on the temperature measurements, observations by a camera and IR camera.

Introduction

Additive manufacturing for building (or printing) 3D structures, e.g., laser cladding process is apparently one of the most promising cost and time intensive technologies. A laser power control is a crucial factor for the high quality of the cladding process (see [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF]). Other factors that influence building a 3D body from a melted powder include the powder supply rate and the velocity of the laser head movements. The process is so complicated that several models of it have been proposed (see, e.g., [START_REF] Picasso | A simple but realistic model for laser cladding[END_REF], [START_REF] Liu | Effects of process variables on laser direct formation of thin wall[END_REF]), but even the simplest of them -adopted from [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF] also in this paper -is described by a nonlinear differential equation, which has to be extended by adding differential equations of a controller and a filter. Thus, the proper selection of the control system structure and tuning of its parameters requires extensive simulations, because direct experiments on a real systems are forbidden as being dangerous -the laser power is about 1 kW. The construction of the robot that performs movements of the laser head is shown in Fig. 1 (left panel), while the right panel shows the laser head with a camera that looks at the cladding process along the laser beam. Notice that in [START_REF] Baraniecki | System for laser microsurfacing of metal powders[END_REF] an original, alternative construction that is specifically designed station, which is based on a three-axis CNC machine equipped with the laser head. The results presented in this paper are expected to be useful also for this new construction.

Our aim in this paper is to propose a methodology for running such simulations. The needs for an extension of control systems simulations methodology stems from the following features of the laser cladding:

the process is a repetitive one, i.e., the laser heating head moves back and forth when developing a 3D object (see [START_REF] Rogers | Control systems theory and applications for linear repetitive processes[END_REF] for a review of repetitive processes and their applications), radiation/optical sensors (e.g. camera, pyrometer) or thermovision -infrared camera can be used for the process monitoring, the repetitive process is time-varying between the passes of the laser head (see the right panel in Fig. 2, where the temperature of a hot, melted metal lake is gradually growing from pass to pass).

Notice that similar features are present in 3D printers and in other additive, powder based technologies, which can make the methodology proposed here more universal.

From the view-point of the designing a control system our in this paper goals are the following:

1. to tune parameters of the PI controller so as to minimize a control quality criterion that is specific for the laser power control task, namely, we penalize too large jumps of the laser power, 2. to decide whether a filter that makes the control signal smoother should be used or not.

The reason for considering the second goal is in extremely large variations in measurements of the lake temperature. This fact can be seen in Fig. 2 (right panel), in which one can observe spikes that differ by ± 50-100 degrees from the mean temperature of the surroundings, while typical errors are at the levels ± 20. These errors are caused either by image processing, since it is very difficult to locate the the hottest place of the lake (see Fig. 2 -left panel) or by a pyrometer that averages the lake temperature, but at the turning points it averages partly the lake temperature and partly the temperature of the surrounding (see [START_REF] Bi | Investigation on the direct laser metallic powder deposition process via temperature measurement[END_REF] for extensive investigations of cladding profess by IR camera). One can expect that such a filter will act also as an anti-wind up module (see, e.g., [START_REF] Mandarapu | Digital pi controller using anti-wind-up mechanism for a speed controlled electric drive system[END_REF] for practical explanations of wind-up phenomenon in control systems with PI controllers). The paper is organized as follows. In the next section we describe the models of:

the dependence of the lake temperature on the laser power, proportional-integral (PI) controller, control signal smoothing module and the methodology of simulations is presented. Then, we use this methodology to reach the two goals mentioned above. Finally, we present an excerpt of the results obtained from experimental implementation of the proposed control system.

Laser Power Control -Model and Simulation Methodology

A mathematical model for one laser pass of the dependence between the laser power Q(t) and the lake temperature y(t) has been adapted from the paper of Tang and Landers Part I [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF] (later abbreviated as TL). It has the following form:

τ y (t) + y(t) = K (Q(t)) β , y(0) = y 0 (1)
where y 0 -given initial temperature. The constants in (1) are defined as follows (compare [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF]):

-β = 6.25 10 -2 is an experimentally selected constant, -τ = 2.96 10 -2 -in sec. is the system time constant.

Overall system gain K (amplification) is defined by

K = K 1 (V α ) (M γ ) = 1413.58, (2)
where:

-K 1 = 0.5 10 β 1.42 10 3 is empirically established system amplification, -V = 3.4 in mm/sec is the laser traverse speed, -M = 4.0 in g/min -the powder supply rate, -α = -7.1 10 -3 and γ = 3.0 * 10 -3 experimentally selected constants.

In [START_REF] Devesse | Design of a model-based controller with temperature feedback for laser cladding[END_REF] a more exact nonlinear model for the lake temperature is derived that takes into account also the laser head position. One can incorporate this model in our simulation methodology, replacing (1) by the model from [START_REF] Devesse | Design of a model-based controller with temperature feedback for laser cladding[END_REF]. However, model (1) occurred to be sufficiently exact for a constant speed of the laser head, since then t simultaneously specifies the head position.

For a multi-pass cladding process we have to extend model (1) by including the heat exchange between passes. We propose the simplest model that can incorporate the influence of the temperature at k-th pass y k (t) on the temperature at (k + 1) pass, i.e., y k+1 (t) of the following form: for passes k = 0, 1 . . .

τ y k+1 (t)+y k+1 (t) = K (W k+1 (t)) β +ξ y k (t), y k+1 (0) = Y k (0), t ∈ (0, T ), (3)
where T > 0 is the pass length (the time that the laser head needs to travel along a 3D object under construction), W k (t) is the laser power at k-th pass at time t, while ξ is the coefficient that governs the influence of the temperature at k-th pass, denoted as Y k (t), on the lake temperature at the next pass. Due to forth and back movements of the laser head, for t ∈ [0, T ] Y k (t) is defined as follows:

Y k (t) = y k (t) if k odd, y k (T -t) if k even. (4)
In order to run (3) we need also the initial condition along the pass, which is assumed to be Y 0

(t) = Y 0 = const, t ∈ [0, T ],
where Y 0 is the temperature of the base.

To check model (3), it was run with a constant laser power W = 0.2 kW. The results are shown in Fig. 3. Comparing it to Fig. 2 (right panel) one can observe similar qualitative behavior, i.e., the lake temperature grows from pass to pass, and then it stabilizes at the level of about 1250-1300 o C. Notice that in Fig. 3 the temperatures from subsequent passes are tilled, while in Fig. 2 they are joined serially.

The analysis of Fig. 2 reveals why a control system for the lake temperature stabilization is necessary. Namely, at the points where the laser head turns back one can observe the essential growth of the lake temperature, since at a short time the laser operates twice at the same regions, near the end points. In [START_REF] Mandarapu | Digital pi controller using anti-wind-up mechanism for a speed controlled electric drive system[END_REF] the proportional-integral (PI) controller is proposed, together by an anti-wind-up by the feedback from the saturation block to the input of PI controller. The control system scheme that is proposed in this paper is depicted in Fig. 4. It consists of PI-controller followed by a control signal smoothing module. The PI-controller that is usually described as follows:

q(t) = K e(t) + 1 T i t 0 e(τ ) dτ , (5)
where K > 0 is the amplification, T i > 0 is the time of doubling of the error signal e(t), while q(t) is the controller output. For the purposes of our simulations it is convenient to rewrite (5) as follows:

q k (t) = K P e k (t) + K I t 0 e k (τ ) dτ, (6)
where K P > 0 is the amplification of P-that, K I > 0 is the amplification of I-path, e k (t) = y ref (t) -y k (t) is the tracking error at k-th pass and y ref (t) is the reference signal, which is assumed to be constant3 in this paper. The integral term in (6) reduces the tracking error and reduces the laser power variability. However, due to large measurement errors, this reduction is not sufficient. For this reason, the control smoothing exponentially weighted moving average (CSEWMA) filter is added in the control loop (see Fig. 4). In fact, CSEWMA is a low pass filter with a specific parametrization that allows its intuitive tuning (see [START_REF] Rafaj Lowicz | Local detection of defects from image sequences[END_REF]). Being the low-pass filter, CSEWMA acts also as an anti-wind-up device by reducing large changes of the laser power. For our purposes the CSEWMA filter is described as follows:

W k (t) = -h W k (t) + h q k (t), W k (0) = W 0 k , (7)
where W 0 k is defined as

W 0 k = W k-1 (0) if k odd, W k-1 (T ) if k even. , (8)
while 0 < h < 1 is a smoothing parameter. Its role is easy to interpret when a finite difference approximation is considered: Hence, the output of the filter at t + ∆ t is a weighted combination of its older output W k (t) and its input q k (t). Smaller values of h force more smoothing.

W k (t + ∆ t) = (1 -h ∆ t) W k (t) + h ∆ t q k (t). ( 9
Algorithm 1 -algorithm for simulating a repetitive control system.

Step 0) Select K P > 0, K I > 0, 0 < h < 1, y ref (t) and the number of passes kmax. Set pass counter k = 1 and

Y 0 (t) = Y 0 , t ∈ [0 T ].
Step 1) Solve the following system of coupled differential and integral equations for t ∈ (0, T )

τ y k (t) + y k (t) = K (W k+1 (t)) β + ξ Y k-1 (t), y k+1 (0) = Y k (0), (10)
q k (t) = K P (y ref (t) -y k (t)) + K I t 0 (y ref (τ ) -y k (τ )) dτ, (11)
W k (t) = -h W k (t) + h q k (t), W k (0) = W 0 k . ( 12
)
Step 2) Update Y k (t) and W 0 k according to, ( 4) and ( 8), respectively. If k < kmax, set k := k + 1 and go to Step 1), otherwise, STOP.

Several remarks are in order concerning the above algorithm.

1. The pass index k is used for y k (t), W k (t) and q k (t) for readability. However, in the implementation it suffices to store only four lists that contain the sampled versions of current W k (t) and q k (t) as well as sampled versions of y new (t) and y old (t) to store the sampled versions of y k (t) and y k-1 (t), respectively, and then to change the contents of y old (t) by y new (t) at Step 2. 2. The algorithm has been extensively tested by the authors. The excerpt of the results is shown in Fig. 6. 3. A general construction of the algorithm is applicable for simulating repetitive control systems for other processes. To this end, it suffices to replace (10) by another one (or by the set of differential equations). 4. The control system in Algorithm 1 is nonlinear. However, the static nonlinearity is rather weak and for verifying the system stability along the pass we can use the results form [START_REF] Rogers | Control systems theory and applications for linear repetitive processes[END_REF], Chapter 7.5 for a linearized system. Our aim in this section is to propose the methodology of running simulations for tuning the control system, i.e., selecting parameters of PI controller and CSEWMA block.

A control quality index J(K I , h, D) depends on parameter K I of PI controller, parameter h of CSEWMA block and on kmax × N matrix D of observation errors. This matrix consists of rows d(

k def = [d (k) 1 , d (k) 2 , . . . , d (k)
N ], where d (k) n 's are realizations of independent random variables with zero mean and a probability distribution selected by the user. For J(K I , h, D) we propose to take the following expression: [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF] where ∆ > 0 is the sampling period (N ∆ = T ), γ > 0 is the weight of the second penalty term, while ŷkmax (t n , d n ) is defined as

J(K I , h, D) = = ∆ N n=1 |y ref (t n ) -ŷkmax (t n , d n )| + γ ∆ N n=2 |W kmax (t n ) -W kmax (t n-1 )|,
ŷkmax (t n , d (kmax) n ) = y kmax (t n ) + d (kmax) n , n = 1, 2, . . . N, (14)
where y kmax (t n ) is obtained from the last pass of Algorithm 1 at n-th sampling instant. Notice that ŷkmax (t n , d

) indirectly depend on all earlier realizations of random errors, since in Algorithm 1 we use the following expression:

qk (t) = K P (y ref (t) -ŷk (t)) + K I t 0 (y ref (τ ) -ŷk (τ )) dτ, (15)
instead of [START_REF] Rafaj Lowicz | Numerical optimal control of integral-algebraic equations using differential evolution with fletchers filter[END_REF], with obvious changes also in [START_REF] Rogers | Control systems theory and applications for linear repetitive processes[END_REF]. In (15) ŷk (t)) is a noise corrupted version of y k (t)) that is obtained as follows: after calculating y k (t) (by solving [START_REF] Rafaj Lowicz | Local detection of defects from image sequences[END_REF]) this signal is sampled and corrupted by random disturbances:

ŷk (t n , d (k) n ) = y k (t n ) + d (k) n , n = 1, 2, . . . , N. (16)
Then, ŷk (t)) is obtained by interpolating ŷk (t n , d

n ). The first term in [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF] is the mean absolute deviation of the system output from the reference value. The second term in ( 13) is less obvious in control theory. It is a penalty for too large changes in a control signal, i.e., the laser power in our case. Notice that this term is similar to, but not the same as, the total variation, which is the norm in the space of functions with bounded variation (see, e.g., [START_REF] Benedetto | Integration and modern analysis[END_REF] for the definition and properties.).

Taking into account that D's is randomly generated, we have to average J(K I , h, D) also with respect to D. The following algorithm calculates the averaged quality index for fixed K I and h

J(K I , h) = j -1 max jmax j=1 J(K I , h, D (j) ), (17)
where D (j) is j-th realization of random matrix D., while j max is the number of samples used for averaging.

Algorithm 2 -algorithm for estimating quality of a repetitive control system for fixed tuning parameters.

Step 0) Select K P > 0, K I > 0, 0 < h < 1, y ref (t) and the number of repetitions j max . Set run counter j = 1 and S = 0.

Step 1) Generate a random matrix D (j) and use its elements in (16).

Step 2) Run Algorithm 1 with (15) replaced by [START_REF] Rafaj Lowicz | Numerical optimal control of integral-algebraic equations using differential evolution with fletchers filter[END_REF].

Step 3) Calculate J(K I , h, D (j) ) according to [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF] and update S = S + J(K I , h, D (j) )/j max .

Step 4) If j < j max , set j := j + 1 and go to Step 1). Otherwise, STOP with S as the output.

Having Algorithm 2 at our disposal we can optimize J(K I , h) with respect to K I and h (and possibly other parameters of the control system). To this end, one can use one of many available and well known global optimization algorithms, e.g., genetic or evolutionary algorithms (see, e.g., [START_REF] Coello | Evolutionary algorithms for solving multi-objective problems[END_REF]) or their recently developed versions that take constraints into account by applying the Fletcher filter [START_REF] Rafaj Lowicz | Fletchers filter methodology as a soft selector in evolutionary algorithms for constrained optimization[END_REF], [START_REF] Rafaj Lowicz | Numerical optimal control of integral-algebraic equations using differential evolution with fletchers filter[END_REF]. The need for using global optimization methods instead of faster but local ones stems from the fact that random errors may lead to false local minima, even if a large number of repetitions j max is used. The time of running Algorithm 2 once, for 6 cladding layers and j max = 600 repetitions, was about 100 sec. on a standard PC i7/3 GHz. In order to grasp an intuition how the landscape of J(K I , h) looks like, Algorithm 2 has been run on the grid for h ∈ [0.1, 0.9] and K I ∈ [0, 10]. The remaining parameters have been selected as follows: K P = K E = 1, ∆ = 0.01 sec., random errors -uniformly distributed in [-20, 20] o C. The results of the simulations are shown in Fig. 5. As one can notice, the landscape is rather flat, but the global optimum is clearly visible. Although a relatively large number of repetitions j max = 600 has been used, one can find many local minima (see, e.g., a vicinity of h = 0.85, K I = 5left panel of Fig. 5). Nevertheless, for the right panel of this figure, one can easily read out that pairs (h, K I ) in the darkest region provide the smallest J(K I , h), which is equal to 4.2. The larger region: h ∈ [0.2, 0.55], K I ∈ [3.5, 8] provides quite good quality of about 4.4. Thus, the cladding process can be run safely in this region.

Selected Experimental Results

The results presented in the previous sections have been based on simulations.

Here, we provide selected experimental results of using the PI-CSEWMA control system to the laser cladding robot.A 3D thin wall has been built several times, using different settings of the control system that are listed in the first and the second column in Table 1. The third column contains J criterion, as defined by [START_REF] Tang | Melt pool temperature control for laser metal deposition processespart i: Online temperature control[END_REF], with λ = 3 and ∆ = 1/6. The reason for selecting this value is the following: the length of each pass was 6 seconds, but for evaluating J we have taken 2000 samples (from 5000 to 7000) in the middle of the pass, which corresponds to one second. The results displayed in the fourth column is the percentage of the improvement of J by introducing control systems in comparison to the case when constant power W = 0.65 kW was applied (first row). As one can observe, introducing the PI controller reduces J by 38 %. Additional reduction, by additional 14 % is obtained when under-smoothed (h = 0.5) CSEWMA filter is introduced. When properly selected smoothing (h = 0.1) in CSEWMA filter is applied, then we obtain the reduction of J to 65 -66 % of its value without control. Notice that interchanging positions of PI and CSEWMA smoother does Table 1. The quality of control (3-rd column) for different control system configurations and tuning (1-st column) and the percentage of the improvement with respect to no control (constant laser power) case (4-th column)

not influence this result. Applying other then zero initial state of CSEWMA filter also does not improve the performance, although it can be desirable in order to get even smoother control signals.

Conclusions and the Directions ff Further Research

The methodology of simulating the dynamics of controlled repetitive processes and optimizing parameters of the control system has been presented in details. Then, this methodology has been applied to the laser power control of the cladding process. The control system consisting of PI-controller and CSEWMA filter has been tuned by the simulations. Then, the designed control system has been tested on a real-life cladding robot. The results are satisfactory in the sense that we have obtained a good performance of the system both in its pass-to-pass behavior as well as in its behavior along each pass, which is smoother then in the case when the CSEWMA filter is not used. Further effort is needed to better control the material deposition at the end points, using images from the IR (thermovision camera), but this topic is outside the scope of this paper.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Cladding robot (left panel) and the laser head with a camera (right panel)

Fig. 3 .

 3 Fig.3. The results of pass-to-pass simulations -the lake temperature growth from pass 1 (solid curve) through pass 2 (the shortest dashed line) to pass 4 (the longest dashed line). The solid straight line -the desired temperature of the lake

Fig. 4 .

 4 Fig. 4. Upper panel -the scheme of laser cladding as a repetitive process. Lower panel -a closed loop system with a controller combined with the CSEWMA module -control along the pass

Fig. 5 .

 5 Fig. 5. Dependence of control quality criterion J(KI , h) on (KI and h -3D plot (left panel) and contour plot (right panel)

Fig. 6 .

 6 Fig.[START_REF] Liu | Effects of process variables on laser direct formation of thin wall[END_REF]. Four passes of the cladding process for approximately optimal tuning of the control system. Left panel -the lake temperature vs time. Right panel -the laser power used

PI+CSEWMA 1

 1 KP = 0.9, KI = 5.88, h = 0.5, q0 = 0 3.95 52PI+CSEWMA 2 KP = 0.9, KI = 5.88, h = 0.1, q0 = 0 2.88 65 PI+CSEWMA 3 KP = 0.9, KI = 5.88, h = 0.1, q0 = 450 2.81 66 CSEWMA + PI KP = 0.9, KI = 5.88, h = 0.1, q0 = 0 2.78 66

It is reasonable to consider also time-varying y ref (t) so as to compensate a shape of 3D body at the laser turning points, but the discussion on this is outside the scope of our paper.


Ryszard Kozera
email: ryszard.kozera@gmail.com

Alexander Prokopenya
email: alexanderprokopenya@sggw.pl

Orthogonal Illuminations in Two Light-Source Photometric Stereo

Keywords: Photometric stereo, image irradiance equation, ambiguity in shape reconstruction, computer vision

In this paper we investigate the case of ambiguous shape reconstruction from two light-source photometric stereo based on illuminating the unknown Lambertian surface. So-far this problem is merely well-understood for two linearly independent light-source directions with one illumination assumed as overhead. As already established, a necessary and sufficient condition to disambiguate the entire shape reconstruction process is controlled by the satisfaction of the corresponding secondorder linear PDE with constant coefficients in two independent variables. This work extends the latter to an arbitrary pair of light-source directions transforming the above constraint into a special nonlinear PDE. In addition, a similar ambiguity analysis is also performed for a special configuration of two light-source directions assumed this time as orthogonal and contained in the vertical plane. Finally, this work is supplemented by illustrative examples exploiting symbolic computation used within a framework of continuous reflectance map model (i.e. an image irradiance equation) and applied to a genuine Lambertian surfaces.

Problem formulation

The task of surface reconstruction from its image (or images) poses one of the major challenges in computer vision. There are two common approaches adopted to tackle this problem. The first one (termed a stereo-like method) involves projective geometry tools applied to multiple camera images [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Faugeras | Three-Dimensional Computer Vision -A Geometric View Point[END_REF][START_REF] Trucco | Introductory Techniques for 3-D Computer Vision[END_REF]. Solving the so--called matching problem constitutes here the gist of the entire reconstruction process based on incorporating the triangulation-like technique [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Faugeras | Three-Dimensional Computer Vision -A Geometric View Point[END_REF][START_REF] Horn | Robot Vision[END_REF]. The second approach, coined shape from shading [START_REF] Trucco | Introductory Techniques for 3-D Computer Vision[END_REF][START_REF] Horn | Robot Vision[END_REF][START_REF] Horn | Shape from Shading[END_REF], relies on a priori knowledge of the physical properties of light reflectance inherent to the specific material coating the unknown surface S. In this model the surface S is assumed to be formed as a graph S = graph(u) of the unknown C k (k = 0, 1) function u : Ω → IR defined over an image domain Ω ⊂ IR 2 . Recall that the normal n(s) = (n 1 (s), n 2 (s), n 3 (s)) to S at the point s = (x, y, u(x, y)) ∈ S (with n 3 (s) < 0) reads as n(s) = (u x (x, y), u y (x, y), -1). Assuming a distant light--source (yielding a parallel light-beam along the direction p = (p 1 , p 2 , p 3 )) the so-called image irradiance equation R(u x (x, y), u y (x, y)) = E(x, y)

encapsulates the relation between reflected light from surface S and absorbed light in the image Ω -see [START_REF] Horn | Robot Vision[END_REF]. Here E(x, y) denotes the intensity of the absorbed light at the image point (x, y) ∈ Ω. On the other hand, R called reflectance map refers to the intensity of the reflected light from s ∈ S with normal n(s). In the special case of a Lambertian surface S L , R L is proportional to cos(α), where α represents the angle between the normal n(s) and a given light-source direction p = (p 1 , p 2 , p 3 ) -see [START_REF] Horn | Robot Vision[END_REF]. Consequently, for S L the corresponding image irradiance equation (1) (over Ω = {(x, y) ∈ IR 2 : 0 ≤ E(x, y) ≤ 1}) reads as: p|n(s) p n(s) = p 1 u x (x, y) + p 2 u y (x, y) -p 3 p 2 1 + p 2 2 + p 2 2 u 2 x (x, y) + u 2 y (x, y) + 1 = E(x, y).

(

) 2
By Cauchy-Schwartz inequality |E(x, y)| ≤ 1 holds. The extra condition required E(x, y) ≥ 0 excludes shadowed (x, y)-points. The equation ( 2) yields the first order non-linear PDE (in two variables) modelling a single image shape from shading for S L . Commonly, in solving (2) one usually searches for u ∈ C k (k = 1, 2), modulo a vertical shift v = u + c (here c is an arbitrary constant). A single image irradiance equation ( 2) renders generically an ill-posed problem (see [START_REF] Brooks | Shading without shape[END_REF][START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF][START_REF] Brooks | Impossible and ambiguous shading patterns[END_REF][START_REF] Bruss | The eikonal equation: some results applicable to computer vision[END_REF][START_REF] Kozera | On complete integrals and uniqueness in shape from shading[END_REF]). Though various mathematical extra assumptions can disambiguate [START_REF] Faugeras | Three-Dimensional Computer Vision -A Geometric View Point[END_REF] or limit its number of solutions [START_REF] Horn | Robot Vision[END_REF][START_REF] Deift | Some remarks on shape-from-shading in computer vision[END_REF][START_REF] Oliensis | Uniqueness in shape from shading[END_REF][START_REF] Chojnacki | A direct computation of shape from shading[END_REF][START_REF] Kimmel | Tracking level sets by level sets: a method of solving the shape from shading problem[END_REF][START_REF] Kozera | Uniqueness in shape from shading revisited[END_REF], they often turn unrealistic or difficult to be met from real camera images. A feasible remedy is to employ a photometric stereo technique for which the unknown surface S L (or S) is consecutively illuminated by distant multiple light--sources positioned along linearly independent directions [START_REF] Horn | Robot Vision[END_REF][START_REF] Woodham | Photometric stereo: a reflectance map technique for determining surface orientation from multiple images[END_REF]. In this set-up, contrary to the classical stereo method relying on images taken from multiple cameras, only a single camera is deployed.

As demonstrated in [START_REF] Horn | Robot Vision[END_REF][START_REF] Woodham | Photometric stereo: a reflectance map technique for determining surface orientation from multiple images[END_REF][START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF], three-light source photometric stereo suffices to guarantee a unique surface S L = graph(u) (up to a vertical shift) with u ∈ C 1 satisfying the corresponding system of three image irradiance equations (2) (formulated for p, q and r) over an image Ω = Ω 1 ∩ Ω 2 ∩ Ω 3 . The respective right-hand sides of (2) (i.e. the intensities E i ≥ 0 over Ω i , for i = 1, 2, 3) are obtained by consecutive illuminations of S L along three linearly independent directions p, q r ∈ IR 3 . The entire reconstruction process is decomposed here into two following steps: a) a gradient computation (an algebraic step):

u x (x, y) = f 1 (x, y, E 1 , E 2 , E 3 ), u y (x, y) = f 2 (x, y, E 1 , E 2 , E 3 ). ( 3
)
In case of three light-sources the resulting vector field ∇u = (u x , u y ) satisfying three equations of type ( 2) is uniquely determined in (3) in terms of image intensities E i (i = 1, 2, 3) and light-source directions p, q and r -see [START_REF] Horn | Robot Vision[END_REF][START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF].

b) a gradient integration (an analytic step):

u(x, y) = u(x 0 , y 0 ) + γ u x (x, y)dx + u y (x, y)dy, (4)
where (x 0 , y 0 ), (x, y) ∈ Ω and γ : [t 0 , t 1 ] → Ω is an arbitrary piecewise-C 1 curve such that γ(t 0 ) = (x 0 , y 0 ) and γ(t 1 ) = (x, y). Note that (x 0 , y 0 ) is fixed here and u(x 0 , y 0 ) represents any constant c. The choice of γ joining (x 0 , y 0 ) with varying (x, y) is arbitrary (at least in continuous model of ( 1)), provided Ω is simply connected and the vector field (3) fulfills the so-called integrability condition.

The latter reads for u ∈ C 2 as

u xy (x, y) = u yx (x, y), (5)
and for u ∈ C 1 as γc u x (x, y)dx + u y (x, y)dy = 0,

holding for any piecewise-C 1 loop γ c ⊂ Ω (with γ c (t 0 ) = γ c (t 1 )). Assuming E i (i = 1, 2, 3) are formed by a genuine Lambertian surface S L = graph(u), a unique vector field ∇u obtained from (3) is automatically integrable and hence the formula (4) determines an unambiguous u ∈ C 1 (up to a constant).

The case of two-light source photometric stereo requires more intricate analysis (see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF][START_REF] Onn | Integrability disambiguates surface recovery in twoimage photometric stereo[END_REF]). The respective system of two image irradiance equations:

p 1 u x (x, y) + p 2 u y (x, y) -p 3 p 2 1 + p 2 2 + p 2 3 u 2 x (x, y) + u y (x, y) 2 + 1 = E 1 (x, y), q 1 u x (x, y) + q 2 u y (x, y) -q 3 q 2 1 + q 2 2 + q 2 3 u 2 x (x, y) + u y (x, y) 2 + 1 = E 2 (x, y) (7)
is solved by the following vector field (see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]):

u x = p (q 1 p|q -p 1 q 2 )E 1 + q (p 1 p|q -q 1 p 2 )E 2 + (p 3 q 2 -p 2 q 3 )ε √ Λ p (p 3 q 2 -q 3 p|q )E 1 + q (q 3 p 2 -p 3 p|q )E 2 + (p 1 q 2 -p 2 q 1 )ε √ Λ , u y = p (q 2 p|q -p 2 q 2 )E 1 + q (p 2 p|q -q 2 p 2 )E 2 + (p 1 q 3 -p 3 q 1 )ε √ Λ p (p 3 q 2 -q 3 p|q )E 1 + q (q 3 p 2 -p 3 p|q )E 2 + (p 1 q 2 -p 2 q 1 )ε √ Λ , (8)
where

Λ = p 2 q 2 1 -E 2 1 (x, y) -E 2 2 (x, y) -p|q ( p|q -2 p q E 1 (x, y)E 2 (x, y)) (9)
with the function ε(x, y) taking values ±1 so that f (x, y) = ε(x, y) Λ(x, y) is continuous (for u ∈ C 1 ) or smooth (for u ∈ C 2 ).

As discussed in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF] an image domain Ω is often decomposed into Ω = Ω (1) ∪ Ω (2) ∪ Γ with Λ > 0 over disjoint sub-domains Ω (j) ⊂ Ω (here j = 1, 2) and with Λ = 0 satisfied along some smooth curve Γ . Consequently, over each component Ω (j) either ε(x, y) ≡ 1 or ε(x, y) ≡ -1 hold, yielding the respective two vector fields (u ±

x , u ± y ) in ( 8) satisfying [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF]. Assuming E 1 and E 2 are generated by a genuine Lambertian surface S L , at least one of these two pairs of vector fields (over each Ω (j) ) is integrable.

The ambiguous case of more than one integrable vector field from ( 8) is well-understood merely for the special configuration of light-source directions, namely when p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ) (here q 2 1 + q 2 2 > 0) -see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]. In fact, generically there is only one integrable vector field over each Ω (j) (and thus over Ω). The latter is governed by the fulfillment of the tight condition [START_REF] Kozera | On complete integrals and uniqueness in shape from shading[END_REF]. Indeed, given two simply connected sub-images Ω (j) , a necessary and sufficient condition enforcing both vector fields (u +

x , u + y ) and (u - x , u - y ) to be integrable (and thus yielding the existence of u + , u -∈ C 2 (Ω (i) ) solving ( 7)) reads as:

q 1 q 2 (u yy (x, y) -u xx (x, y)) + (q 2 1 -q 2 2 )u xy (x, y) = 0, (10)
which is to be satisfied by either u + or u -. In addition (see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]), if one of u + (or of u -) satisfies ( 10) then so does the complementary one i.e. u -(or u + ).

In a rare situation of ( 10) holding, both pairs (u + , u -), determined by (4) over Ω (j) (j = 1, 2), can bifurcate (i.e. can be glued together) along the curve Γ to yield either zero or two or four global solutions u ∈ C 2 over the whole image Ω (see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]). As it turns out, there are extra geometrical relations between both graphs of u + and u -. Indeed, let K u (x, y) denote the Gaussian curvature of the surface S L = graph(u) taken at the point (x, y, u(x, y)) ∈ S L and determined by the formula (see [START_REF] Do Carmo | Differential Geometry of Curves and Surfaces[END_REF]):

K u (x, y) = u xx (x, y)u yy (x, y) -u 2 xy (x, y) (1 + u 2 x (x, y) + u 2 y (x, y)) 2 . (11)
Interestingly, it is proved in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF] that if both ( 7) and ( 10) hold then:

K u + (x, y) = -K u -(x, y). (12)
Consequently, if s 1 = (x, y, u ± (x, y)) ∈ S ± L yields a hyperbolic point (having negative Gaussian curvature) then s 2 = (x, y, u ∓ (x, y)) ∈ S ∓ L renders an elliptic point (having positive Gaussian curvature). Thus for p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ) the convexity/concavity ambiguity is automatically excluded. Noticeably, such ambiguity eventuates for single image shape from shading with p = (0, 0, -1) (see e.g. [START_REF] Horn | Robot Vision[END_REF][START_REF] Brooks | Impossible and ambiguous shading patterns[END_REF]). Finally, if s 1 is a parabolic point (with vanishing Gaussian curvature) then so is s 2 .

This paper extends the above results with the following (see Section 2):

1. The necessary and sufficient condition [START_REF] Kozera | On complete integrals and uniqueness in shape from shading[END_REF] for testing the ambiguity in ( 7) is extended to arbitrary pairs of linearly independent light-source directions p = (p 1 , p 2 , p 3 ) and q = (q 1 , q 2 , q 3 ) Ultimately, such general case leads to the non-linear second-order PDE (see [START_REF] Kozera | On shape recovery from two shading patterns[END_REF]), which inherently constitutes a challenging problem.

2. To alleviate the above difficulty, special configurations of two light-source directions (different from p = (0, 0, -1) and q = (q 1 , q 2 , q 3 )) are here admitted. Namely, the family of orthogonal unit vectors (parallel to OXZ-plane):

p α = (-1 -α 2 , 0, -α), q α = (α, 0, -1 -α 2 ) (13)
(with 0 < α < 1) is considered. The corresponding analysis for solving [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF] (with p = p α and q = q α ) addresses the uniqueness issue and establishes, in rare ambiguity cases, intrinsic geometrical inter-relations between multiple solutions in two light-source photometric stereo (see also Th. 1).

The theoretical Section 2 is also supplemented (see Section 3) by illustrative examples supported by Mathematica numerical and symbolic computation applied to the continuous Lambertian model with pixels represented as ideal points (x, y) ∈ Ω and image intensities (E 1 , E 2 ) simulated here from u as left-hand sides in [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF]. Finally, the closing Section 4 summarizes the main thrust of this paper together with indicating its extensions and hints their possible solutions.

Orthogonal Illumination Directions in Vertical Plane

In the first part of this section (i.e in (i)) we consider an arbitrary configuration of linearly independent light-source directions p and q with p = q = 1. Under such general assumption, condition (10) (testing whether the second vector field from (8) over each Ω (j) is also integrable) is subsequently extended to the nonlinear PDE expressed by [START_REF] Kozera | On shape recovery from two shading patterns[END_REF]. Noticeably, the latter constitutes a difficult task for further theoretical analysis. In order to deal with the latter somehow, a special case of orthogonal vectors p and q contained in the OXZ-plane is here admitted (see (ii)). The analysis to follow complements already established results in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF], covering a different special choice of p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ) with q 2 1 + q 2 2 > 0. (i) Assume now that E 1 and E 2 introduced in (7) are generated by a genuine Lambertian surface S L = graph(u), where u ∈ C 2 . Evidently, still by [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF] both image intensities are also expressible in terms of u x and u y . Using symbolic computation in Mathematica [START_REF]Wolfram Mathematica[END_REF][START_REF] Budzko | Symbolic-numerical methods for searching equilibrium states in a restricted four-body problem[END_REF][START_REF] Budzko | Symbolic-numeric analysis of the equilibrium solutions in the restricted four-body problem[END_REF] (or alternatively see complicated proof in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF]) one arrives at:

Λ(x, y) = (u x (p 2 q 3 -p 3 q 2 ) -u y (p 1 q 3 -p 3 q 1 ) + p 2 q 1 -p 1 q 2 ) 2 u 2 x + u 2 y + 1 = ( p × q|n ) 2 u 2 x + u 2 y + 1 . ( 14
)
Combining the latter with (8) (over each Ω (j) determined by Λ > 0) leads to:

u x = u x , if sgn(ε)sgn(θ) > 0; (a 2 -b 2 -c 2 )ux+2acuy+2ab 2abux+2bcuy+b 2 -a 2 -c 2 , if sgn(ε)sgn(θ) < 0 (15)
and

u y = u y , if sgn(ε)sgn(θ) > 0; 2acux+(c 2 -a 2 -b 2 )uy+2bc 2abux+2bcuy+b 2 -a 2 -c 2 , if sgn(ε)sgn(θ) < 0, (16)
where

a = p 3 q 2 -p 2 q 3 , b = p 1 q 2 -p 2 q 1 , and c = p 1 q 3 -p 3 q 1 , (17)
and θ = u x (p 2 q 3 -p 3 q 2 ) -u y (p 1 q 3 -p 3 q 1 ) + p 2 q 1 -p 1 q 2 . Here the function ε = ε(x, y) is everywhere constant taking values ±1. The first pair (u x , u y ) from ( 15) and ( 16) satisfies integrability condition (5) over each Ω (j) as u ∈ C 2 . On the other hand, the integrability of the second vector field from ( 15) and ( 16) eventuates, if and only if:

(a 2 -b 2 -c 2 )u x + 2acu y + 2ab 2abu x + 2bcu y + b 2 -a 2 -c 2 y = 2acu x + (c 2 -a 2 -b 2 )u y + 2bc 2abu x + 2bcu y + b 2 -a 2 -c 2 x .
Upon resorting to the symbolic computation in Mathematica, the last equation is transformable into the following non-linear PDE:

c (a -bu x (x, y)) u yy (x, y) + a 2 -c 2 + bcu y (x, y) -abu x (x, y) u xy (x, y) +a (bu y (x, y) -c) u xx (x, y) = 0, (18)
which is to be satisfied by u. Evidently, the latter does not hold generically and therefore the resulting integrability condition [START_REF] Kozera | On shape recovery from two shading patterns[END_REF] disambiguates almost always two-source photometric stereo modelled by [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF] -see also Ex. 1. Still, in a pursue of solving a rare ambiguity in [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF], one ought to deal with [START_REF] Kozera | On shape recovery from two shading patterns[END_REF] which inevitably leads to a non-trivial task. Thus in the next step (ii) of this section, a tighter constraint imposed on illumination directions p and q is considered.

(ii) Suppose now that two light-source directions p α and q α are introduced according to [START_REF] Chojnacki | A direct computation of shape from shading[END_REF]. The resulting two image irradiance equations coincide with:

α - √ 1 -α 2 u x (x, y) u 2 x (x, y) + u y (x, y) 2 + 1 = E 1 (x, y), αu x (x, y) + √ 1 -α 2 u 2 x (x, y) + u y (x, y) 2 + 1 = E 2 (x, y). (19)
Since p α ⊥ q α (i.e. are orthogonal), the function Λ defined in (9) simplifies into:

Λ(x, y) = 1 -E 2 1 (x, y) -E 2 2 (x, y) = u 2 y (x, y) 1 + u 2 x (x, y) + u 2 y (x, y) . ( 20
)
Thus

Ω Λ>0 = {(x, y) ∈ Ω : u y (x, y) = 0} and Ω Λ≡0 = {(x, y) ∈ Ω : u y (x, y) = 0}.
Similarly to the special case of p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ) discussed in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF], often the triples u, p α and q α (see Ex. 1) yield Ω Λ>0 = Ω (1) ∪ Ω (2) , with Ω (j) (j = 1, 2) standing for two disjoint sub-domains of Ω and Ω Λ≡0 = Γ representing a smooth curve in Ω. Here a bifurcation curve (along which solutions over Ω (1) and Ω (2) are glued) coincides with the curve (overlapping with Ω Λ≡0 ):

Γ = {(x, y) ∈ Ω : u y (x, y) = 0}. ( 21
)
Furthermore in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF] , since p 2 = q 2 = 0 then a = b = 0 and since q 1 = -p 3 = α and p 1 = q 3 = -√ 1 -α 2 then c = 1. Consequently, both formulae (15) and ( 16) are reducible into:

u x =      u x , if sgn(ε)sgn(θ) > 0; u x , if sgn(ε)sgn(θ) < 0; u y =      u y , if sgn(ε)sgn(θ) > 0;
-u y , if sgn(ε)sgn(θ) < 0. ( 22) Thus for the non-generic case of the second vector field (u x , -u y ) also integrable (over each Ω (i) ), the function u should satisfy the following linear PDE:

u xy (x, y) = 0. ( 23
)
Note that the last equation can be independently reached by substituting a = b = 0 and c = 1 into [START_REF] Kozera | On shape recovery from two shading patterns[END_REF]. The generic case of ( 23) not fulfilled is illustrated also in Ex. 1. Upon double integration (first over x and then over y) of ( 23), the following representation for u (if two vector fields from [START_REF] Budzko | Symbolic-numerical methods for searching equilibrium states in a restricted four-body problem[END_REF] are to be integrable over

Ω (i) ) holds: u(x, y) = φ 1 (x) + ψ 1 (y), (24)
for some twice continuously differentiable functions φ 1 and ψ 1 in a single variable. It is not difficult to show that φ 1 (x) = u(x, 0) -c, where ψ 1 (0) = c. Similarly, ψ(y) = u(0, y) -φ 1 (0) = u(0, y) -u(0, 0) + c. Naturally, an analogous argument applies to the second solution v ∈ C 2 to (19) resulting in v satisfying [START_REF] Budzko | Symbolic-numeric analysis of the equilibrium solutions in the restricted four-body problem[END_REF] and hence v(x, y) = φ 2 (x) + ψ 2 (y), where φ 2 and ψ 2 are defined similarly to the introduction of φ 1 and ψ 1 . Furthermore, by [START_REF] Budzko | Symbolic-numerical methods for searching equilibrium states in a restricted four-body problem[END_REF], the second function v ∈ C 2 (over each Ω (j) ) fulfills (v x , v y ) = (u x , -u y ). Combining the latter with u x (x, y) = φ 1 (x), u y (x, y) = ψ 1 (y), v x (x, y) = φ 2 (x) and v y (x, y) = ψ 2 (y) yields "a conjugate-like" relation between u and v:

u(x, y) = φ 1 (x) + ψ 1 (y), v(x, y) = φ 1 (x) -ψ 1 (y) + c 1 , (25)
with c 1 being a constant. Formula (25) determines specific analytic representations of the solutions to [START_REF] Onn | Integrability disambiguates surface recovery in twoimage photometric stereo[END_REF] over each Ω (j) in a rare situation of the ambiguous two light-source photometric stereo. A straightforward verification shows that the geometrical constraint from ( 12) is also preserved for p α and q α . Indeed, combining [START_REF] Deift | Some remarks on shape-from-shading in computer vision[END_REF] with [START_REF] De Boor | A Practical Guide to Spline[END_REF] y). Note that the condition (21) coupled with [START_REF] De Boor | A Practical Guide to Spline[END_REF] implies that Γ (or more general Ω Λ≡0 ) represents a line L = {(x, y * ) ∈ Ω : ψ 1 (y * ) = 0} (or a collection of lines) parallel to the OX-axis. In addition, as ∇u = (u x , u y ) and ∇v = (u x , -u y ) visibly any critical point of u is also a critical point of v.

leads to K u (x, y) = (φ 1 (x)ψ 1 (y))(1 + (φ 1 (x)) 2 + (ψ 1 (y)) 2 ) -2 = -K v (x,
Moreover, if such point represents a local minimum (maximum, saddle) for u then it is also a local maximum (minimum, saddle) for v. Note also that by [START_REF]Wolfram Mathematica[END_REF] any critical point of u (and thus of v) belongs to the set Ω Λ≡0 and thus to the potential bifurcation curve Γ . The non-generic ambiguity case discussed above for p α and q α is illustrated in Ex. 2. Evidently, as [START_REF] Kozera | On shape recovery from two shading patterns[END_REF] or [START_REF] Budzko | Symbolic-numeric analysis of the equilibrium solutions in the restricted four-body problem[END_REF] are generically not fulfilled, there exists only one solution u ∈ C 2 to (7) (or [START_REF] Onn | Integrability disambiguates surface recovery in twoimage photometric stereo[END_REF]) over each Ω (j) . Upon gluing u together, only one global solution u ∈ C 2 prevails over entire image Ω. On the other hand, the rare scenario of the existence of two solutions u, v ∈ C 2 over each Ω (j) (i.e. satisfying [START_REF] Budzko | Symbolic-numeric analysis of the equilibrium solutions in the restricted four-body problem[END_REF]) leads to the possible bifurcations along Γ rendering 0, 2 or 4 global solutions of class C 2 over entire image Ω. Recall that we assume here Ω = Ω (1) ∪ Ω (2) ∪ Γ . The other decomposition topologies of Ω = Ω Λ>0 ∪ Ω Λ=0 are discussed in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF]. The detailed analysis justifying necessary and sufficient conditions to guarantee successful C k bifurcations (k = 0, 1, 2) (in case of Ω = Ω (1) ∪ Ω (2) ∪ Γ ) holding along Γ exceeds the scope of this paper and therefore is here omitted.

Taking into account the above argument, the main theoretical contribution of this paper can be summarized into the following: Theorem 1. Assume that u ∈ C 2 (which graph S L = graph(u) represents an illuminated genuine Lambertian surface) satisfies [START_REF] Onn | Integrability disambiguates surface recovery in twoimage photometric stereo[END_REF] (or (7)) with p α and q α determined by (13) (or with arbitrary linearly independent p and q). Suppose, moreover that Λ in (9) satisfies Λ > 0 over simply-connected Ω (j) (j = 1, 2). In order that there exists only one more solution v ∈ C 2 to (19) over Ω (j) it is necessary and sufficient for u to satisfy (23) (or [START_REF] Kozera | On shape recovery from two shading patterns[END_REF]). In addition, the ambiguous case for p α and q α yields u(x, y) = φ(x) + ψ(y) and v(x, y) = φ(x) -ψ(y), with φ and ψ determined as in [START_REF] Kozera | Curve modeling via interpolation based on multidimensional reduced data[END_REF]. Finally, the Gaussian curvatures of graphs of u and v at respective points (x, y, u(x, y)) and (x, y, v(x, y)) (see [START_REF] Deift | Some remarks on shape-from-shading in computer vision[END_REF]

) satisfy K u (x, y) = -K v (x, y).
As already pointed out, a local ambiguity to [START_REF] Onn | Integrability disambiguates surface recovery in twoimage photometric stereo[END_REF] or [START_REF] Brooks | Circularly symmetrical eikonal equations and non-uniqueness in computer vision[END_REF] (over Ω Λ>0 ) can even be more proliferated to a global one (over entire image Ω) due to possible bifurcations of u and v along Ω Λ≡0 . The matter gets more complicated if Ω Λ≡0 forms an open subset of Ω (see [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]). The latter occurs once p × q|n = 0 as implied by [START_REF] Kimmel | Tracking level sets by level sets: a method of solving the shape from shading problem[END_REF]. Due to these intricacies, the respective discussion on bifurcation issue in two light-source photometric stereo is here left out.

Experiments

This section includes two examples illustrating the main results established in Section 2 (see also Th. 1). The experiments presented here are carried out with the aid of Mathematica symbolic computation. The corresponding pictures of images (with the respective intensities E 1 and E 1 ) are simulated synthetically upon admitting arbitrary or specific illumination directions and assuming a genuine u as temporarily initially given.

Example 1. (a) Consider a Lambertian hemi-sphere

S 1 L = graph(u 1 ) with u 1 ∈ C 2 ( Ω) defined as u 1 (x, y) = R 2 -x 2 -y 2 , where Ω = {(x, y) ∈ IR 2 : x 2 +y 2 ≤ R 2 }.
For two linearly independent normalized light-source directions p and q the respective image irradiance equations read as:

p 1 u x + p 2 u y -p 3 u 2 x + u 2 y + 1 = E 11 (x, y) = -p 1 x -p 2 y -p 3 R 2 -x 2 -y 2 R , q 1 u x + q 2 u y -q 3 u 2 x + u 2 y + 1 = E 21 (x, y) = -q 1 x -q 2 y -q 3 R 2 -x 2 -y 2 R , (26)
over

Ω = Ω 1 ∩ Ω 2 ∩ Ω, where Ω i = {(x, y) ∈ IR 2 : E 1 i1 (x, y) ≥ 0} (for i = 1, 2)
. The negative values of E i1 represent shadowed subareas of Ω. The bifurcation curve Γ (see ( 14)) reduces into a planar quadratic determined by:

(p 3 q 2 -p 2 q 3 )x -(p 3 q 1 -p 1 q 3 )y + (p 2 q 1 -p 1 q 2 ) R 2 -x 2 -y 2 = 0. (27)
Furthermore, the condition (18) ascertaining the existence of exactly one solution u 1 ∈ C 2 (Ω) (modulo its vertical shift) to ( 26) enforces u 1 to satisfy:

(cx -ay) ax + cy -b R 2 -x 2 -y 2 = 0, (28)
over Ω (j) , with a, b and c defined as in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF]. Clearly, the equation ( 28) is not satisfied by u 1 for all (x, y) ∈ Ω (j) . This yields uniqueness in solving (26) (i.e. u = u 1 ) within u ∈ C 2 (Ω). In particular, for p α and q α (see [START_REF] Chojnacki | A direct computation of shape from shading[END_REF]), ( 17) combined with a = b = 0 and c = 1 reduce ( 27) into y = 0 rendering Γ = {(x, y) ∈ Ω : y = 0}, Ω (1) = {(x, y) ∈ Ω : y < 0} and Ω (2) = {(x, y) ∈ Ω : y > 0}. The ambiguity condition ( 28) is transformed into xy = 0 merely fulfilled along both X-and Y -axes. Hence, again for arbitrary p α and q α , uniqueness of u 1 prevails. This is expected since p α and q α represents a special case of general positions of p and q analyzed above.

(b) Let a Lambertian hill-like surface S 2 L = graph(u 2 ) with u 2 ∈ C 2 ( Ω) be defined according to u 2 (x, y) = (2(1 + x 2 + y 2 )) -1 , over e.g. Ω = {(x, y) ∈ IR 2 : |x| ≤ 1 and |y| ≤ 1}. The respective two image irradiance equations read as:

p 1 u x + p 2 u y -p 3 u 2 x + u 2 y + 1 = E 12 (x, y) = -p 1 x -p 2 y -p 3 (x 2 + y 2 + 1) 2
x 2 + y 2 + (x 2 + y 2 + 1) 4 , q 1 u x + q 2 u y -q 3 u 2

x + u 2 y + 1 = E 22 (x, y) = -q 1 x -q 2 y -q 3 (x 2 + y 2 + 1) 2

x 2 + y 2 + (x 2 + y 2 + 1) 4 , (29)
over unshadowed Ω = Ω 1 ∩Ω 2 ∩ Ω, where Ω i = {(x, y) ∈ IR 2 : E 1 i2 (x, y) ≥ 0} (for i = 1, 2). The bifurcation curve Γ from ( 14) is defined by the following equation:

(p 2 q 1 -p 1 q 2 ) + (p 3 q 2 -p 2 q 3 )x 2(1 + x 2 + y 2 ) 3/2 + (p 3 q 1 -p 1 q 3 )y 2(1 + x 2 + y 2 ) 3/2 = 0. ( 30
)
On the other hand, the integrability condition (18) stipulates u 2 to satisfy

(cx -ay) b + 6(ax + cy) 1 + x 2 + y 2 = 0, (31)
over Ω (j) , with a, b and c introduced as in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF]. Again, (31) is not fulfilled by u 2 for all (x, y) ∈ Ω (j) . Hence there exists a unique solution of class C 2 to (29) (i.e. u = u 2 ) over Ω. In the special case of p α and q α (see [START_REF] Chojnacki | A direct computation of shape from shading[END_REF]), [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF] coupled with a = b = 0 and c = 1 transform (30) into y = 0. Hence again Γ = {(x, y) ∈ Ω : y = 0}, Ω (1) = {(x, y) ∈ Ω : y < 0} and Ω (2) = {(x, y) ∈ Ω : y > 0}. Furthermore the ambiguity constraint (31) is reduced into xy = 0 which holds again only along both X-and Y -axes. Thus for p α and q α the uniqueness of u 2 to (29) eventuates (as also follows from the above general (p, q)-position analysis).

Example 2. Consider now a Lambertian paraboloid S 3 L = graph(u 3 ) with u 3 ∈ C 2 ( Ω) defined as u 3 (x, y) = (x 2 +y 2 )/2 over Ω = {(x, y) ∈ IR 2 : |x| ≤ 1 and |y| ≤ 1}. For a general position of p and q, the integrability condition [START_REF] Kozera | On shape recovery from two shading patterns[END_REF] reads here as b(ay -cx) = 0. Provided b = 0, such constraint never holds over any open subset of Ω and thus uniqueness of u 3 follows. Noticeably the bifurcation curve Γ (a line) is determined here by p 2 q 1 -p 1 q 2 + (p 2 q 3 -p 3 q 2 )x + (p 3 q 1 -p 1 q 3 )y = 0.

For p α and q α from ( 13) the resulting two image irradiance equations are:

α - √ 1 -α 2 u x (x, y) u 2 x (x, y) + u y (x, y) 2 + 1 = E 13 (x, y) = α - √ 1 -α 2 x x 2 + y 2 + 1 , αu x (x, y) + √ 1 -α 2 u 2 x (x, y) + u y (x, y) 2 + 1 = E 23 (x, y) = αx + √ 1 -α 2 x 2 + y 2 + 1 , (32)
over Ω α = Ω (for α = 1/ √ 2) and over Ω α = {(x, y) ∈ Ωα : -1 ≤ x ≤ 0.204} (for α = 1/5). Fig. 1 shows images of S 3 21) coincides here with the X--axis (i.e. here y = 0) trimmed either to -1 ≤ x ≤ 1 or to -1 ≤ x ≤ 0.204, for α = 1/ √ 2 or α = 1/5, respectively. The remaining decomposition components of Ω α (along which Λ > 0) read as Ω (1) 2. Furthermore, since by [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF] the constant b = 0, the ambiguity condition b(ay -cx) = 0 is now satisfied. Hence, upon combining [START_REF] De Boor | A Practical Guide to Spline[END_REF] with φ(x) = x 2 /2 and ψ(y) = y 2 /2, the only one another C 2 solution to (32) equals to u 4 (x, y) = (x 2 -y 2 )/2 (modulo a vertical shift), over each Ω (j) α (here j = 1, 2). Once u k are glued with itself (k = 3, 4) along X-axis, two C 2 class global solutions to (32) over entire image Ω α are defined (i.e. u 3 and u 4 ) -see Fig. 3 for α = 1/ √ 2 (with gluing curve). On the other hand the local solutions u 3 (or u 4 ) over Ω α to yield next two only C 1 class solutions u 34 and u 43 to (32) over Ω α . Indeed, a C 2 differentiability is excluded as lim (x,y)→(x,0) u 3yy (x, y) = 2 = -2 = lim (x,y)→(x,0) u 4yy (x, y). The remaining two C 1 global solutions to (32) over Ω 1/ √ 2 (i.e. u 34 and u 43 ) are plotted in Fig. 4. The case for α = 1/5 differs merely by different Ω 1/5 . Note that by [START_REF] Deift | Some remarks on shape-from-shading in computer vision[END_REF] the respective Gaussian curvatures yield K u3 (x, y) = 4(1 + x 2 + y 2 ) -2 = -K u4 (x, y) which is consistent with Th. 1. A unique critical point (0, 0) of u 3 , u 4 , u 34 and u 43 belongs, as expected to Γ .

L with α = 1/ √ 2 over Ω 1/ √ 2 = Ω = [-1, 1] × [-1, 1]. The bifurcation curve Γ α in (
α = {(x, y) ∈ Ω α : y > 0} and Ω (2) α = {(x, y) ∈ Ω α : y < 0} -see Fig.

Conclusion

This paper extends the claims of [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF], where a special configuration of two light-source directions p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ) (for q 2 1 + q 2 2 > 0) in photometric stereo is studied. In particular, a respective integrability condition expressed by the non-linear PDE [START_REF] Kozera | On shape recovery from two shading patterns[END_REF] is derived here for a general configuration of linearly independent p and q. Not unexpectedly, such PDE forms a difficult theoretical problem and consequently a new family of orthogonal illumination directions p α and q α is introduced in (13). The corresponding system of two image irradiance equation ( 19) is subsequently analyzed and the resulting ambiguity versus uniqueness question is addressed in Th. 1. As proved in this paper, similar ambiguity results are obtainable for p α and q α defined in (13) (see Th. 1 and Section 1) as compared to those already established in [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF] for p = (0, 0, -1) and q = (q 1 , q 2 , q 3 ). The experiments reported in Section 3 illustrate the main results from Section 2. A possible extension of this paper includes a relevant ambiguity analysis performed for the general positions of two linearly independent light-source directions. Inevitably, any pending argument would rely on characterizing the multiple solutions to the non-linear PDE defined in [START_REF] Kozera | On shape recovery from two shading patterns[END_REF]. Another worth investigation venue is to complete a global uniqueness-ambiguity analysis for two image photometric stereo over entire image Ω (see also [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Kozera | On shape recovery from two shading patterns[END_REF]). The latter should first cover the case of p α and q α defined in [START_REF] Chojnacki | A direct computation of shape from shading[END_REF] and then should refer to an arbitrary selection of two linearly independent illumination directions. In particular, an extra attention should be paid here in derivation of analytical and numerical methods designed to localize possible bifurcation curve(s) Γ . Some relevant clues concerning this task can be found in [START_REF] Kozera | Curve modeling via interpolation based on multidimensional reduced data[END_REF][START_REF] De Boor | A Practical Guide to Spline[END_REF][START_REF] Kozera | External versus internal parameterization for length of curves with nonuniform samplings[END_REF][START_REF] Kvasov | Methods of Shape-Preserving Spline Approximation[END_REF][START_REF] Kozera | Piecewise-quadratics and exponential parameterization for reduced data[END_REF][START_REF] Atamanyuk | The algorithm for optimal polynomial extrapolation of random processes[END_REF]. Finally, a similar analysis based on real Lambertian images, where pixelization and noise occur, forms another vital extension topic of this paper. Usually, handling noisy and digitized image data requires a more robust integration techniques dealing computationally with u recovery phase determined by theoretical formula (4). Such analytic step is often substituted by the pertinent optimization task designed to compute numerically the closest discrete integrable vector field -for more see e.g. [START_REF] Horn | Robot Vision[END_REF][START_REF] Horn | Shape from Shading[END_REF][START_REF] Noakes | Non-linearities and noise reduction in 3-source photometric stereo[END_REF].

( 1 )

 1 α cross-bifurcate along Γ 1/ √ 2 with u 4 (or with u 3 ) over Ω[START_REF] Faugeras | Three-Dimensional Computer Vision -A Geometric View Point[END_REF]

Fig. 1 .√ 2 and q 1

 121 Fig. 1. Two images of S 3 L = graph(u3) illuminated along p 1/ √ 2 and q 1/ √ 2 directions.

Fig. 2 .α∪

 2 Fig. 2. Decomposition of Ωα = Ω (1) α ∪ Ω (2) α ∪ Γα for (32) with a) α = 1/ √ 2, b) α = 1/5.

Fig. 3 .√ 2 =Fig. 4 .

 324 Fig. 3. Two C 2 class global solutions u3 and u4 to (32), over Ω 1/ √ 2 = [-1, 1] × [-1, 1].