Advanced Image Processing Techniques for Physics Studies T. Craciunescu and A. Murari

with contribution from:

G. Kocsis, P. Lang, I. Tiseanu, J. Vega and JET EFDA Contributors

*

*See the Appendix of F. Romanelli et al., Fusion Energy Conference 2008 (Proc. 22nd Int. FEC Geneva, 2008) IAEA, (2008)

Workshop on Fusion Data Processing Validation and Analysis, ENEA- Frascati (26-28 March 2012)

Optical flow - extraction of advanced information for control and physics studies

- Implementation (CLG, MPEG) - Application to pellets and instability tracking

Automatic instability detection

- Phase congruency image classification - Sparse image representation for disruption prediction - Interest points and local features for image identification

Optical flow

attempt to find the vector field which describes how the image is changing with time frame t optical flow frame ( t+1) Basic assumptions: the grey values of image objects in subsequent frames do not change over time

f

(

x



u

,

y



v

,

t

 1 ) 

f

(

x

,

y

,

t

)  0 small displacements:

f x u



f y v



f t

 0 Ill-posed problem: small perturbations in the signal can create large fluctuations in its derivatives undetermined set of equations

Combined local-global (CLG) method

Assumes that the unknown optic flow vector is constant within some neighbourhood of size ρ.

 A sufficiently large value for ρ is very successful in rendering the method robust under noise.

 in flat regions where the image gradient vanishes, the become again undetermined.

problem 

Coarse-to-fine multi-resolution approach E CLG

     1 

w T J

   3

f



w

    2   

dxdy

Incorporates a global smoothness assumption for the estimated flow field.

 Larger values for α result in a stronger penalisation of large flow gradients and lead to smoother flow fields.

 At locations with | ∇ f| ≈ 0, no reliable local flow estimate is possible, but the regulariser | ∇ u|2 + | ∇ v|2 fills in information from the neighbourhood the filling-in effect.

Predictive understanding of the underlying processes of the pellet-plasma interaction

Recent investigations revealed that pellet ablation is a complex 3D process taking place on the μs timescale → pellet cloud dynamics (expansion, instabilities and drifts) Analysis of pellet cloud dynamics and drifts by observing the visible radiation with fast framing cameras algorithms and by applying image processing __________________________________________________________________________________

* detailed results will be presented at EPS conference (G. Kocsis et al) Experiments with sophisticated diagnostic settings performed during the 2011 campaign of AUG

Determination of ice extrusion velocity by optical flow method

Image sequences provided by a CCD camera viewing the ice at the exit of the nozzles of the extrusion cryostat Illustration of optical flow calculations showing the extruded deuterium ice in case of JET pulse #76379 Line profiles through the images and its reconstruction (bottom)

MPEG-2 compressed space

• Statistical redundancies in both temporal and spatial directions: - inter-pixel correlation - simple translation motion between consecutive frames

(I)- Intrinsic frames

- coded using only information present in the picture itself by d iscrete cosine transform ( DCT) – – Processing at the level of MB8 blocks DCT concentrates the energy into the low frequency coefficients (spatial redundancy) → → neglecting the low value coefficients High-frequency coefficients are more coarsely quantized than the low-frequency coefficients (P)- Predicted frames are coded with forward motion compensation, using the nearest previous reference (of type I or P) images.

→ Motion is represented by a field of motion vectors (MV) → one MV per macroblock (B) - Bi-directional frames are also motion compensated, this time with respect to both past and future reference frames.

 The parts of the image that do not change significantly are simply copied from other areas or other frames.

 In case of the other parts, for each MB16, the best matching block is searched in the reference frame(s).

Encoding macroblocks is implemented using MB14

MV field used:   as a crude initial estimation for optical flow recovery for image segmentation  Confidence measure to ensure that the MV field is meaningful Assumption: areas with strong edges exhibit better correlation with real motion than textureless ones  weighted averages of the image gradients can be expressed using DCT coefficients: • • eigenvalue decomposition size of the eigenvalue is a measure of uncertainty in the direction of the corresponding eigenvector (the stronger the eigenvalue, the lower the error variance)

Error estimation

 Peak signal-to-noise ratio (PSNR) residual image > 14 dB of the  The difference between the speeds of the different pellets in the same ribbon structure below 12.5%.

Computing time

Image processing step Time (ms) Segmentation using information from MPEG video compressed domain

Applying the regularization rules for the 0.6

MV field Median filtering + 0.8

Morphological operations

Optical flow calculation performed on the segmented image region

Image derivatives SOR iteration Median filtering 1.8

4.1

0.3

 Total optical flow computation time:

~16.4 ms

.

 Image acquisition framing rate:

50 Hz

 Optical flow evaluation can be engrafted in MPEG compressing routines in case of real time estimation of the speed of moving objects.

Tracking of plasma instabilities

MARFEs can reduce confinement leading to harmful disruption → a risk for the integrity of the devices MARFEs determine a significant increase in impurity radiation → a clear signature in the video data

Image processing step

Applying the regularization rules for the MV field Median filtering Dilation/Erosion Labeling Objection centroid determination

Time (ms)

0.5

0.7

3.0

2.4

0.3

Phase congruency

Visually discernable features coincide with those points where the Fourier waves, at different frequencies, have congruent phases Extraction of highly informative features at points of high PC

Mach bands

construction of PC from the Fourier components

black

– measued luminance

red

–brightnesses as perceived

Lateral inhibition vs.

Pahse congruency

M.C. Morrone et al., Mach bands are phase dependent, Nature 324(1986)250.

t Approximations of

F

and

G

by convolving the signal with a quadrature pair of filters (linear-phase filters for phase information preservation) 

Q n even

,

n Q odd

  symmetric/antisymmetric quadrature pairs of nonorthogonal wavelets An appropriate choice for constructing the symmetric/antisymmetric quadrature pairs of filters are the Gabor filters.

Gabor filters with different frequencies and orientations Response for Gabor filter oriented vertically

Combine all the orientations SIM map, pooled into a single similarity score  96.2% were correctly interpreted. From the misclassified events 0.03% were false positives and 3.5% false negatives.

Sparse learned representations of video images

Sparse image representation

Since images are usually large, the decomposition is implemented on overlapping patches instead of whole images.

 D – fixed, general (DCT, wavelet) or it can be adapted to suit the application domain.  Learning both D and  in an efficient way has been the focus of much of recent published work.

• • D initialized from random patches of natural images.

Then learned adaptively from the data such that the decomposition is sparse:

D

, min

i i

 1 ,...,

k



i x i

 

i D

2 2

s

.

t

.

 0 

L

Each patch written as a column vector Sparsity: - counting the number of non-zero elements in a vector • • •

Matching pursuit algorithm

first find the one D atom that has the biggest inner product with the signal then subtract the contribution due to that atom repeat the process until the signal is satisfactorily decomposed

 Decomposition error:

R



X

,

D

 

X

 

D

2   The dictionary trained on patches from natural images 

N

different classes

S i

of signals – learn separate dictionaries, one per class – a signal belonging to one class is reconstructed poorly by a dictionary corresponding to another class.

– classification is performed by using residual reconstruction errors of a signal by the dictionary belonging to a class

R



x

,

D i

 as a discriminative operator for classification.

Limited results: ~ 85 % success classification rate

Learning discriminative dictionaries

‘Good’ for one class ‘bad’ for the orher by incorporating discriminative components:

D

, min

i i

 1 ,...,

k



i



x i

 

i D

2 2   

i

1   

i

 

j D i T D j

2

F A F

 

i j a ij

2 → → atoms representing common features in all classes tend to appear repeated almost exactly in dictionaries corresponding to different classes False similar reconstruction decomposition error • • Detect by inspecting the inner product of dictionary atoms.

Threshold for controlling the sharing atoms Dictionary incoherence term

Encourages dictionaries associated to different classes to be as independent as possible, while still allowing for different classes to share features.

Improved (preliminary) results: ~ 92 % success clasification rate   Further tuning adjustable parameters – mainly size of patches multiscale Multiscale framework to capture first a global appearance of objects

Image retrieval

 Bag-of-words model - represention of a ‘sentence’ as an unordered collection of words, disregarding grammar and even word order.

 Transforms an image into a large collection of feature vectors invariant to: • • • image translation scaling rotation • • illumination changes local geometric distortion → – algorithms to detect and describe local features in images: MSER (Efficient Maximally Stable Extremal Region)  extremal region

R i



p



R i

, 

q



boundary

 

i

 Im  Im  

Component tree

 Rooted, constructed connected by tree successive thresholdings account taking hierarchic into image inclusion → maximally stable regions are those regions which have approximately the same region size across 2Δ neighboring threshold images Features calculated: mean gray value, region size, center of mass, width, dimension of the bounding box (weights for the features can be used to adapt to different kinds of input data).

Matching criteria: smallest Euclidean distance between feature vectors

Reference image Image identification using maximally stable extremal regions Various MARFE images

Conclusions

 Optical flow method for the study of several fusion plasma relevant issues, able to provide the real velocity for objects moving close to structures.

– – MPEG Motion segmentation - a key contrivance to allow very fast optical flow estimation Application to pellet injection and pellet dynamic  Phase congruency as a highly localized operator for automatic MARFE identification with a good prediction rate.

 Sparse image representation for disruption prediction. Encouraging preliminary results. Improvements expected mainly from a more efficient definition and implementation of the discriminative operator.

 Image retrieval by image local feature detection.