On the Quick Convergence of Arbitrarily

Shaped PSF Estimation for Single Image

Blind Debluring

Submitted By
Engr. Muhammad Waleed

Supervisor
Dr. Aftab Khan

DEPARTMENT OF COMPUTER SYSTEMS ENGINEERING (DCSE)

UNIVERSITY OF ENGINEERING AND TECHNOLOGY (UET)
PESHAWAR, PAKISTAN

August2017
iii
 On the Quick Convergence of Arbitrarily
Shaped PSF Estimation for Single Image
Blind Debluring

By

Engr. Muhammad Waleed

A Thesis
presented to the University of Engineering and Technology (UET), Peshawar in partial
fulfilment for the degree requirement of

Master of Sciences

in

Computer Systems Engineering

2017

Department of Computer Systems Engineering (DCSE)

University of Engineering and Technology (UET)
Peshawar, Pakistan
iv
 DECLARATION

I hereby declare that I am the sole author of this thesis. This is a true copy of the thesis,
including any required final revisions, as accepted by my examiners. It is further
declared, that I have fulfilled all the requirements in line with the Quality Assurance
guidelines of the Higher Education Commission.

I understand that my thesis may be made electronically available to the public.

v
 ABSTRACT
Over the past few decades’ large technological advances have taken place, Nowadays, most of
image data are captured from handheld and mobile phone cameras, satellites, CCTVs and other
imaging sources. Images often suffer from degradation due to imperfection in the capturing and
imaging process; with a recorded image inevitably emerging as a degraded version of the
original image. These degradations are caused from various sources like lens defocus, optical
imperfections in the case of a digital camera or atmospheric blurring in the case of satellite/aerial
photography, etc. In real life when we captured images on the camera, after capturing the images
the original scenes are no more but we have left just with that capture information. Now in this
specific case no information about the impulse response of the images is available and in turn
leads to the restoration of these degraded images blindly and is termed as Blind Deconvolution.
To recover such images many algorithms has been was developed like Total Variant,
Richardson-Lucy, Wiener filters etc. Image Quality Measures(IQMs) has been also developed
for low quality images, which is not refer for the high-quality images.

The research does not deal with the parametric PSF and denoising (noise removal) of the image
may it be artificial or deblurring noise. Also, ringing effects of an image is out of scope for the
research. Undoing these imperfections, in some cases, is crucial for many image processing and
multidimensional signal processing tasks. Digital Image Restoration is the different operations
applying on images to achieve the clean approximation of actual image. In many fields like
Defence, Medical, Planetary Sciences, Satellite Imaging etc., the key thing is to obtain the clean
informative image from a blur image. The aim of this research was to investigate into methods
which can aid in the faster convergence of the blind image restoration scheme.

This research used different heuristic optimization algorithms for single image deblurring to
solve the problem. The main contribution of this research is a study of three heuristic algorithms
Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO) and Genetic Algorithm
(GA) and compare empirically the convergence of these schemes. The experimental result and
performance comparison of the algorithms shows that PSO has much better efficiency than ACO
and GA. The implementation of PSO and ACO algorithm is complicated as compared to GA but
the results are efficient.

6
 ACKNOWLEDGEMENTS
First of all, thanks to Almighty Allah (Subhanahu Wa-ta’ala) Who blessed me with the gift of
knowledge, determination and gave me the strength to carry out this research work in such a
pleasant way. I wish to express deep gratitude to my supervisor Dr. Aftab Khan at the
Department of Computer Systems Engineering (DCSE), University of Engineering and
Technology(UET) Peshawar, for his supervision, guidance and valuabletime polishing me
towards higher studies. He taught me like a strict teacher but alwaystreated me like a best friend.
More than a teacher, he will always be a best friend to me.

I want to give a special thanks to Dr. Nasir Ahmad for always helping and guiding me whenever
I had problems in my studies. He has always lent a helping hand greeted with a humble smile. I
am greatly thankful to the close support and guidance of Dr. Laiq Hasan, Chairman DCSE, UET
Peshawar. As the department chair, my teacher and mentor, he stressed on the importance of
higher studies and developing interpersonal skills.

I am thankful to my brother Muhammad Dawood Khan who made my time at home special and
aided a lot in making research work and thesis write-up enjoyable. I love his great company and
caring nature. I dedicate my education to the skillfulness of my teachers, love, care and
encouragement of my family and friends and their devotion of precious time in advising me.

7
Dedicated to my loving parents for their endless love and care
and for their selfless guidance, encouragement and support

8
Table of Contents
Abstract ......................................................................................................................................................... 6
Acknowledgements....................................................................................................................................... 7
List of Publications ..................................................................................................................................... 11
List of Figures ............................................................................................................................................. 12
List of Tables .............................................................................................................................................. 13
List of Acronyms ........................................................................................................................................ 14
1. INTRODUCTION ................................................................................................................................... 16
1.1. Background ..................................................................................................................................... 16
1.2. Research Scope ............................................................................................................................... 20
1.3. Aims................................................................................................................................................. 20
1.4. Motivations ..................................................................................................................................... 20
1.5. Objectives........................................................................................................................................ 22
1.6. Research Methodology ................................................................................................................... 22
1.7. Thesis Overview .............................................................................................................................. 23
2. Image Deconvolution And Deblurring ................................................................................................ 25
2.1. Introduction .................................................................................................................................... 25
2.2. Research Problem ........................................................................................................................... 25
2.3. Blurring Kernels ............................................................................................................................... 26
2.4. Gaussian Blur .................................................................................................................................. 27
2.5. Motion Blur ..................................................................................................................................... 29
2.6. Out-of-focus Blur............................................................................................................................. 30
2.7. Restoration Techniques .................................................................................................................. 32
2.7.1. Inverse Filter ............................................................................................................................... 32
2.7.2. Least-Squares Filters ................................................................................................................... 33
2.7.3. Iterative Blind Deconvolution (IBD) Technique .......................................................................... 35
2.7.4. Richardson Lucy Algorithm ......................................................................................................... 36
2.7.5. Regularization Method ............................................................................................................... 37
2.7.6. Total Variation Algorithm............................................................................................................ 38
3. Overview of Optimization algorithms ................................................................................................. 39
4. Optimization Schemes ........................................................................................................................ 45
4.1. Genetic Algorithm (GA) ................................................................................................................... 45
4.2. Ant Colony Optimization (ACO) ...................................................................................................... 46

9
4.3. Particle Swarm Optimization (PSO) ................................................................................................ 50
5. EXPERIMENTAL RESULTS FOR OPTIMIZATION OF PSF ESTIMATION .................................................. 52
5.1. Experimental Setup and Results ..................................................................................................... 52
6. DISCUSSION AND ANALYSIS ................................................................................................................ 57
7. CONCLUSIONS AND FUTURE WORK ................................................................................................... 58
Conclusions ............................................................................................................................................. 58
Research Findings ................................................................................................................................... 58
Future Work ............................................................................................................................................ 59
References ................................................................................................................................................... 60
Appendices.................................................................................................................................................. 65
A. Appendix A ..................................................................................................................................... 65

10
LIST OF PUBLICATIONS

• Muhammad Waleed, Aftab Khan and Ashfaq Khan, “On the Quick Convergence of
PSF Estimation for Single Image Blind Deblurring”, IEEE International
Conference on Electronics, Computers and Artificial Intelligence, 29 June - 01 July
2017.

11
LIST OF FIGURES
Figure 1.1: Real life blurred images.examples .............................................. Error! Bookmark not defined.
Figure 1.2: Research methodology ................................................................ Error! Bookmark not defined.
Figure 2.1: Blurring Model ......................................................................................................................... 26
Figure 2.2:Model of image degradation and restoration ............................................................................. 28
Figure 2.3:(a) Perspective plot of a Gaussian PSF of size 15x15 with variance σ=2.5 and (b) its respective
frequency domain representation ................................................................................................................ 29
Figure 2.4: (a) Blurred video frame (b) Deblurred using estimated atmospheric turbulence PSF.............. 29
Figure 2.5: (a) Perspective plot for a PSF of linear motion blur of length 15 and angle 0 degrees and (b)
its Fourier transform. .................................................................................................................................. 30
Figure 2.6: A frame from the video of a bar-code depicting vertical motion blurring. This particular frame
was extracted afterwards. ............................................................................................................................ 30
Figure 2.7: (a) Perspective plot of camera out of focus PSF for R=17, and (b) its Fourier transform........ 31
Figure 2.8: (a) Original book image and (b) Out-of-focus blurred version captured using Microsoft
LifeCam. ..................................................................................................................................................... 32
Figure 2.9: Block diagram of Iterative Blind Deconvolution (IBD) algorithm [5,6].................................. 36
Figure 3.1:Input image and restored images with IPSO-DWT and PSO-DWT. .......... Error! Bookmark not
defined.
Figure 3.2:(a) Original Image (b) Restored Image with PSO-DWT (c) Restored Image with IPSO-DWT
....................................................................................................................... Error! Bookmark not defined.
Figure 3.3:Comparison of PSO-DWT and IPSO-DWT................................. Error! Bookmark not defined.
Figure 3.4: Images with different deblurring techniques (a) Original Image (b) Blind deconvolution (c)
Wiener filter (d) Lucy-Richardson (e) Regularized filter (f) Proposed method ........... Error! Bookmark not
defined.
Figure 3.5: Different deblurring techniques and Improvements .................... Error! Bookmark not defined.
Figure 3.6: Graph shows the CPU time for optimal solution. Continues solid line shows absence of the
algorithm to the problem. Dot line shows the algorithm optimal solution .... Error! Bookmark not defined.
Figure 4.1:GA stages of execution................................................................. Error! Bookmark not defined.
Figure 4.2:ACO stages of execution .............................................................. Error! Bookmark not defined.
Figure 4.3: (a) Ants find out the path between nest and food source. (b) when obstacle comes in their
path: then Ants decide by their intelligent behavior whether to turn left or right. (c) More ants pass
through the shortest path so the pheromone amount is stronger than other paths. (d) Ant follow the path
with much amount of pheromone and by this way they choose the shortest path [1]... Error! Bookmark not
defined.
Figure 4.4: PSO stages of execution .............................................................. Error! Bookmark not defined.
Figure 5.1: (a)Real motion blurred image(b) Deblurred using GA (c) Deblurred using ACO (d)Deblurred
using PSO scheme.......................................................................................... Error! Bookmark not defined.
Figure 5.2: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO ....................................................................... Error! Bookmark not defined.
Figure 5.3: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO ....................................................................... Error! Bookmark not defined.
Figure 5.4: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO ....................................................................... Error! Bookmark not defined.

12
LIST OF TABLES
Table 5.1: Results of ACO, GA and PSO ...................................................... Error! Bookmark not defined.

13
List of Acronyms

IQMs Image Quality Measures

TV Total Variation
RL Richardson Lucy
PSF Point Spread Function
DWT Discrete Wavelet Transform
BID Blind Image Deblurring/Deconvolution
PSNR Peak Signal to Noise Ratio
IPSO Improved Particle Swarm Optimization
Eqn. Equation
Fig. Figure
GUI Graphical User Interface
HVS Human Visual System
IBD Iterative Blind Deconvolution
HPF High Pass Filtering
SNR Signal to Noise Ratio
NASA National Aeronautics and Space Administration
JAXA Japan Aerospace Exploration Agency
ESA European Space Agency
CNSA China National Space Administration
IR Impulse Response
TF Transfer Function
DIR Digital Image Restoration
CT Computed Tomography
MRI Magnetic Resonance Imaging
VI Visual Information
CCD Charge Coupled Device
WF Wiener Filter

14
LPF Low Pass Filter
GF Gaussian Filter
OTF Optical Transfer Function
COC Circle of Confusion
CLSF Constrained Least-Squares Filter
LSC Least-Squares Filter
IBD Iterative Blind Deconvolution
FFT Fast Fourier Transform
IFFT Inverse Fast Fourier Transform

15
1. INTRODUCTION
1.1. BACKGROUND

Digital image restoration is a broad field of engineering for advance level research to recover the
original image from a degraded one. Various digital image restoration applications and
techniques have been developed to overcome this challenging problem. Digital image restoration
aims to achieve an approximation of the pristine image which was lost due to noise, miss-focus,
camera-handshake etc.[2, 3].

In medical science, molecular spectroscopy, planetary sciences, satellite imaging, biomedical

images [4, 5], aerial and satellite photography [6-8], handheld camera pictures and videos [2, 9-
11], audio and seismic signals [12, 13], remote sensing data [14], industrial tomography
images[15], astronomical photos[16, 17]and other sources of signals and images the key thing is
to approximate the original image from a blur and noisy image. Over the past few decades, large
technological advances have taken place in the field of image processing with images being
captured by handheld and mobile phone cameras, satellites, CCTVs and other imaging sources.

Technological advancements have taken place all over the world in different fields specially in
engineering. Like in digital image processing area images has been taken for advance level
research, as different agencies with advance technologies in the world are carrying out various
researches at space like National Aeronautics and Space Administration (NASA), Japan
Aerospace Exploration Agency (JAXA), European Space Agency (ESA), China National Space
Administration (CNSA)[18] and many more using imaging data for their research. Most of
blurred image data is taken from handheld and mobile phone cameras, satellites, CCTVs and
other imaging sources. Digital Image Restoration is a broad field of engineering for advance
level research to recover the original image from a degraded one.

Images often suffer from degradation due to imperfection in the capturing and imaging process;
with a recorded image inevitably emerging as a degraded version of the original image. These
degradations are caused from various sources like lens defocus, optical imperfections in the case
of a digital camera or atmospheric blurring in the case of satellite/aerial photography, etc.

16
Undoing these imperfections, in some cases, is crucial for many image processing and
multidimensional signal processing tasks. Lot of digital image restoration applications comes in
front of us now a day. Digital Image Restoration is actually the different operations applying on
images to get the specific information or to achieve the clean approximation of actual image
which is lost due to miss-focus or motion blur (Blur which creates due to motion).

Until now, the PSF for uniform blurring has been estimated using the functional form for some
common types of blur i.e. the Gaussian blur, motion blur and out-of-focus blur. In real life,
especially in the case of motion blur resulting from camera handshake, the blur follows
convoluted paths resulting in complex shaped PSFs. There exist many schemes in the literature
dealing with the restoration of such images with some listed in [2, 3, 9, 19]. Their shape cannot
be easily modelled by a simple equation or defined by a mathematical model for a set of its
parameter(s).

PSFs in such cases exhibit an arbitrary shape and deblurring in this case requires intricate ways
to estimate the blur and effectively remove it. Arbitrarily shaped PSFs have been shown to exist
in the case of atmospheric turbulence blur as well [20]. The focus of the study to investigate into
methods which can help in the quick convergence of the blind image restoration scheme
visualized for arbitrarily shaped PSF estimation scheme.

In restoration of blind images the actual problem is that’s of lack the information about the
images. In daily life when we want to capture the image, it’s all about the specific scene or object
or specific information that we want to capture. For example, the bird sitting on the tree or a
butterfly with a beautiful flower, after getting or capturing these scenes most probably the scene
will be no more like you may say that the bird fly and no more on the tree or the butterfly flew
away and no more with the flower so the only information we have is that the captured data.

But as we already discussed the sometime the image get blur due to motion, environmental noise
etc. then in that case we have to restored the original image which is quite challenging due to the
lack of information like impulse response or transfer function and that’s what we call it blind
image deblurring.

If the captured image is represented by g, it can be modelled as the convolved version of the
original image f and the impulse response of the capturing medium represented by h.

17
 g = f *h + n (1.1)

Here n represents the noise in the system which is assumed absent in our research work. The
recovery or approximation of the original image requires the deconvolution of the blurred image
by using a deblurring filter such as Wiener filter, Richardson-Lucy filter, Total Variation (TV)
filter etc.[21-23]. Furthermore, the recovery requires the approximation of the blurring function
or the Point Spread Function (PSF).

This is challenging as the original scenery information is not available on the image is captured
and the capturing medium’s PSF is not easy to approximate [24]. PSF estimation techniques are
usually slow as the search space is large and they are further marred by the presence of large
image data especially multi-channel images.

To reduce the time taken to quickly approximate the true PSF, parallelization of the algorithm
can be done but this is not encouraging. Rather, optimization techniques are sought which aim at
converging quickly and estimating the PSF fast. In daily life, blurred images come in front of us
taken by an immature photographer have no knowledge of photography or handshake or miss
focus. Big question for the researcher to remove all these distortions and to apply inverse
scenarios to achieve the right goal are all related to the field of Advance Image Restoration.
Fig.1.1 show the real life blurred images.

18
Fig.1.1: Real life blurred images Examples[20]

19
Advance research work in Digital Image Restoration was started in late 1960s while scientist was
studying the space. During the study, different images obtained was disfigured. disfiguring in the
images were caused due to spinning, vibration, trembling and movement of the space craft. Like
space in many other field the meaningful information was distorted in the imaging data and it
becomes difficult to achieve the goal due to loss of information in the imaging data. Here comes
the time which compel research community to work on advance level on imaging
data[25].Besides from the astronomical images recovering the field of image restoration is now
alsowidely using in computer vision, medical imaging and remote sensing [4, 26] etc.

The degradation model or blur kernel of the images in image restoration is known and it is called
Point Spread Function (PSF). When no information given about the image as the original scene
is no more and just we have the degraded image with lack of information. In that scenario, we
restored the image blindly and such process is termed as “Blind Image Deconvolution (BID)”.
It’s a technique for estimating scene and PSF of the image [17].

1.2. RESEARCH SCOPE

The research investigated blind deblurring of space invariant blur for single image restoration.
Space invariant blur includes both parametric and non-parametric (arbitrarily) point spread
functions (PSF). In space-invariant blur we will deal with the arbitrarily shaped PSF.

The research does not deal with the parametric PSF and denoising (noise removal) of the image
may it be artificial or deblurring noise. Also, ringing effects of an image is out of scope for the
research.

1.3. AIMS

The aim of this research is to investigate into methods which can help in the quick convergence
of the blind image restoration scheme visualized for arbitrarily shaped PSF estimation scheme.

1.4. MOTIVATIONS

Motivational factors for the blind image deconvolution in the field of digital image restoration
are as follow:

20
• Medical Imaging: lot of work has been done in digital image processing in medical
sciences.Medical sciencesarehighly facilitated through advance technology like using
imaging data.In Medical, imaging data is using to explore the inner human body like X-Ray
images for different parts of the human bones, ultrasound images also called sonography
which is uses high frequency sound waves, Computed tomography (C.T) scan uses x rays to
take the images if different parts of the body, MRI (Magnetic Resonance Imaging) images,
mammographic images and many others commonly used for diagnosing different type of
human body damages like orthopedic damage (related to bones of the human body), foreign
objects, and tumors etc. in most cases the restoration for all these images is necessary to
achieve clear results.

• Computer imaging: Computer imaging is an exited area now a day related to the visual
information. In modern age of communication, the visual information is transmitted in the
form of digital image. This is the most important area in current race of technology. As we
know the primary sense is a visual sense and when the information processed by the
computer and conveyed through image is remember for the centuries. A lot of information
in words could be kept in a single image. As just one picture says thousands of words. Now
a day as massive data need to be transfer so the field of computer imaging further extend to
various sub fields as image compression and segmentation etc. the two main categories of
the computer imaging field are

▪ Computer vision

▪ Image processing

• Astronomical Imaging: Digital image restoration area also broadly using for the
astronomical imaging. Astronomical images are normally taken with electronic detectors
like CCD (Charge Coupled Device). Detectors like this also found in digital cameras.

As discussed, the image restoration idea is using everywhere when restoration of images
necessary for the information. In most cases the astronomical images effected by noise and
blur and thus the original information vanishes or not observable. Blind image
deconvolution is better for recovering the original astronomical images from the degraded

21
 images as in this case it is not possible to get the prior information about the original image
scene.

• Video Surveillance: In most cases the video streams that we got from surveillance
cameras are blurred or noisy and one can’t figure out the required information from that
stream easily. Therefore, it is required to clean up these video streams from blurring and
noise so, that we get our required data/information from the stream.
These are not the complete applications of the digital image processing.

1.5. OBJECTIVES

The main objectives of the research

• To develop a scheme which quickly converge on true PSF for arbitrarily shaped PSF
estimation
• To investigate evolutionary algorithms apart from genetic algorithm for search
optimization
• To investigate the convergence time of various optimization schemes

1.6. RESEARCH METHODOLOGY

The proposed methodology is

• Intensive literature review

• Analyze various optimization schemes
• Implement the optimization schemes
• Compare with benchmark results
• Validate results
• Disseminate research work through international conference/journal

22
 To find and Study different
Optimization Schemes

To Implement the Schemes

Simulation for Comparison

Simulation in Matlab

Analysis for the best Results

Fig.1.2 Methodology

1.7. THESIS OVERVIEW

Chapter 2 describes the imaging formation processes, some commonly blurring PSFs (causing
the distortion of an image), recovery processes of these degraded images (Image Restoration),
blind deconvolution and Image Quality Measures (IQMs). Chapter 2 has discussed some of the
well-known BID schemes and Blind IQMs in details. Thisinformationis then used as a base to
tackle the image restoration problem.

23
Chapter 3 presented a novel approach towards the BID problems. In this approach image fusion
technique is utilized for the reducing the effort of the tuning for the optimum parameters values
for the restoration filters. The proposed scheme “Fused Restoration Filtering for Blind Image
Deblurring” successfully validated for artificial as well as real life blurred images.

Chapter 4 explained different non-reference IQMs for gauging the quality of the restored
images with different BID schemes. Different blind IQMs such as Spatial Kurtosis, BRISQUE,
NIQE, SSEQ and CQA are discussed and used as quality assessment for restored images.

Chapter 5 discusses and analyses the results for blind IQMs and BID schemes employed in this
research work. It also provides the experimental results of the new proposed algorithm for BID.

Chapter 6summarizes and concludes the research and gives direction for further research. It also
discusses the major contributions of this research study towards the field of BID.

Appendices present the images’ details used in this research, widespread outcomes of different
techniques for the images used in these simulations of the deblurring assessment over a large set
of images.

24
2. IMAGE DECONVOLUTION AND DEBLURRING
2.1. INTRODUCTION

This chapter will discuss the modeling of the image formation processes and some of the
restoration methods. Some of the mathematical and signal processing schemes are proposed for
the BID problems. Also, the overview of the basic concept of the bling deblurring has been
presented in this chapter. This includes the problems, different type of blurs and restoration
filters which are specifically used for the restoration techniques.

2.2. RESEARCH PROBLEM

Digital Image Restoration technique comprise some form of relation between the blurring kernel;
that make a connection between the original scene before the image taken and the degraded
image after capturing the scene. The resultant image is the convolution of the original scene
convolve with the PSF along with some noise. If the captured image is represented by g, it can be
modelled as the convolved version of the original image f and the impulse response of the
capturing medium represented by h.

g = f *h + n (2.1)

Here n represents the noise in the system which is assumed absent in this research work. The
recovery or approximation of the original image requires the deconvolution of the blurred image
by using a deblurring filter such as Wiener filter, Richardson-Lucy filter, Total Variation (TV)
filter etc.

Fig.2.1: Blurring Model

25
To write the above equation Eq.2.1 in frequency domain as followed:

G = FH + N (2.2)

This research work does not deal with the denoising (noise removal) of the image may it be
artificial or deblurring noise. Also, ringing effects of an image is out of scope for the research.

Here both the deblurring and inherent noise are assumed to be absent. In case of noiseless images
the above equation Eq.2.2 become as:

G
F = or F  = GH −1 (2.3)
H

Here H −1 is the inverse filter and F'is the estimate of the original image after restoration.

Model of image degradation and restoration process is shown as:

f
Original
Image
g
Blurred Restoration
Image filter f'
h
Impulse
Response or Noise Restoration Process
PSF
n
Degradation Process

Fig.0.1: Model of image degradation and restoration

As we discussed already g is the degraded result of the original and impulse response (PSF)
along with noise n which is not the concern of this research. f’ is the restored image through
restoration filter.

2.3. BLURRING KERNELS

Blurring can be found in any imaging process which utilizes electromagnetic radiation like X-
rays and visible light. Diffraction restricts the quality of an imaging device to features on the
order of the illuminating wavelength. Scattering of light between the target object and imaging

26
system (like, by the changing refractive index of atmosphere) induces some extra blurring.
Lenses and mirrors cause blurring because they have limited spatial extent and optical faultiness.
Discretization results in yet more blurring because devices such as Charge Couple Devices
(CCDs) average lighting over regions rather than sampling it at discrete points. The blurring
functions are usually related to the following two classes:

• Space invariant

In this class of PSFs the common form of blur PSFs are independent of image pixel
location. These PSF functions produce an identical blurring effect for each pixel
location.

• Space variant

In this class of PSFsthe blurring PSFs, thatcreates an unalike blurring effect depending
on image pixel location. This results in the blurring effect being different for different
pixels.

PSFs can be divided into two types based on shape, namely parametric and non-parametric:
• Parametric form
This form of PSF can easily be modeled as a function and can be expressed by a
mathematical form or equation.
• Non-parametric form
These non-parametric PSFs have arbitrary shapes, therefore these are also called
arbitrary shaped PSF and they cannot be easily modeled in mathematical form. Images
degraded with such type of PSFs are a challenging task.
Blurred images with three different type of blurring will be restored in this research work. These
are discussed in full details in the following sections.

2.4. GAUSSIAN BLUR

The Gaussian blur effect is produced in an image when it is filtered by a low pass filter estimated
by a 2-D Gaussian function. It is also called Gaussian smoothing function. The Gaussian filter in
two dimensions over PSF’s rows and columns, m and n, according to [27] is given as

27
 2 2
−( m +n )
1
h ( m, n ) = e 2
2
(2.4)
2 
2

According to the [28, 29], the PSF for atmospheric turbulence blur can be described by Eq. 2.4,
which can be viewed as a Gaussian blurring, Fig. 2.3 shows the Gaussian blur PSF of size 15x15
and its respective Optical Transfer Function (OTF). The OTF or Fourier transform
approximation of the Gaussian PSF is also a Gaussian function. Fig. 2.4(a) shows the effect of
atmospheric turbulence blur (approximated by Gaussian blur). The low pass filter blurs the
image especially the edges (high pass signal). The image when restored by deblurring becomes
crisper with the edges recovered.

(a) (b)
Fig. 0.2: (a) Perspective plot of a Gaussian PSF of size 15x15 with variance σ=2.5 and (b) its
respective frequency domain representation.

(a)

28
 (b)
Fig. 0.3: (a) Blurred video frame (b) Deblurred using estimated atmospheric turbulence PSF.

2.5. MOTION BLUR

Motion blur is the superficial freckling of speedily moving objects in an image or in a movie or
animation. It is caused by relative motion between the imaging device and the recorded scene.
This relative motion can be translational, rotational, an abrupt change in the scale, or a
combination of all or some of them. In this research work, only the universal translational motion
is considered.When the scene to be recorded translates relative to the camera at a constant
velocity under an angle of φradians with the horizontal axis during the exposure interval, the
distortion is one-dimensional. Denoting the length of motion by L, the angle by φ, the PSF is
given by Eq. 2.5 with reference to [28].

1 L m
 if m2 + n2  and = − tan 
h(m, n; L,  ) =  L 2 n (2.5)
0 elsewhere

Herem and n are the PSF pixel’s coordinates. Fig. 2.5(a) shows the PSF obtained with
application of Eq. 2.5 for linear motion for length of 15 pixels and at an angle of zero degrees
while its spectral domain representation is shown in Fig. 2.5(b). Fig. 2.6 shows the effect of
motion blurring video frame of the bar-code image. The filter spreads the effect of the
neighboring pixels in the direction of motion.

29
 (a) (b)
Fig.0.4: (a) Perspective plot for a PSF of linear motion blur of length 15 and angle 0 degrees and (b)
its Fourier transform.

Fig.0.5: A frame from the video of a bar-code depicting vertical motion blurring. This particular
frame was extracted afterwards.

2.6. OUT-OF-FOCUS BLUR

When a camera takes a picture of scene or body, some parts of the scene are not in the proper. If
the aperture of the camera is round, then the image of any point source is a small disk, referred as
the Circle-Of-Confusion (COC). The measure of the defocus (the diameter of the COC) is
dependent on the focal length and the aperture value of the lens, and the distance between the
object, body or scene and the image capturing device. A perfect model describes the diameter of

30
the COC as well as intensity distribution within the COC accurately. However, if the measure of
the defocusing is large relative to the wavelengths considered, a geometrical approach can be
followed resulting in a uniform intensity distribution within the COC. The spatially continuous
out-of-focus blur of radius R, with PSF coordinates m and n, is given by Eq. 2.6 with reference
to[28].
 1
 if m 2 + n 2  R 2
h(m, n; R) =  CR 2
(2.6)
0 elsewhere

HereC is a constant that must be chosen so that the energy conservation law is satisfied. Fig. 2.7
(a) shows the original PSF and Fig. 2.7 (b) shows its spectral domain representation. One can
notice the low pass behavior (in this case both horizontally and vertically) in Fig. 2.8, as well as
characteristic pattern of spectral zeros in Fig. 2.7 (b).

(a) (b)

Fig.0.6: (a) Perspective plot of camera out of focus PSF for R=17, and (b) its Fourier transform.

Fig. 2.8 shows an image of the MATLAB book and it’s out of focus version captured by
manually changing the focus away from the focal point on the Microsoft LifeCam. The small
size text becomes unreadable as a result of defocusing.

31
 (a)

(b)

Fig.0.7: (a) Original book image and (b) Out-of-focus blurred version captured using Microsoft
LifeCam.

2.7. RESTORATION TECHNIQUES

Different type of restoration techniques is discussed here in this section which are analyzed and
utilized in experiments in this research work.

2.7.1. INVERSE FILTER

If we look onto the Eq. 2.3, it can be easily deduced that the deblurring problem is an ill-posed inverse
problem. So, ideally to undo the blurring the inverse of the PSF is required. If PSF is known than
degraded images can easily be restored. If there is no noise in the blurred image than the deblurring
can be achieved easily in the frequency domain [5]. In frequency domain, the convolution converted
into multiplication. The inverse filtering process can be expressed as;

32
 N
F = F + (2.7)
H

However, in most cases especially real blurred scenarios the PSF are not available.
Unfortunately, many problems exist with inverse filtering, i.e. may inverse filter not exist
because H is zero at selected frequencies. Secondly, if the blurring function’s spectral form
Hdoes not go actually to zero but becomes small, than the second term in Eq. 2.7, (known is
inverse filtered noise) will become very large, thus inverse filtered images are therefore often
dominated with by excessively amplified noise (ill-conditionedness/ill-posedness)[28].

2.7.2. LEAST-SQUARES FILTERS

To address the problem of noise sensitivity/implications associated to the inverse filtering, a

number of restoration techniques are invented and these are collectively called Least-Squares
Filters. In this research work two mostly used methods of Least-Squares Filters are discussed in
the following sections.

• Weiner Filter
Weiner filter is the linear spatially invariant filter, in which the PSF is chosen so that it
minimizes the Minimum Square Error (MSE) between the original and restored images.
It uses statistical approach towards image deblurring; it removes both additive noise and
blurring simultaneously from the blurred images. In frequency domain, it can be
expressed as;

H*
F ' wen = (2.8)
Sn
H +
2

Sf

HereH* shows complex conjugate of H, Sfrepresents power spectral of original image

andSnshows power spectral of the noise.

Let

K= S n
(2.9)
S f

33
Then Eq. 2.8 can be written as;

H*
F ' wen = (2.10)
H +K
2

This K term is also known as the Noise to Signal Ration (NSR). The excessive noise
amplification is not present here. If we put K= 0 than the above Eq. 2.10 becomes the
inverse filter. However, residual blur and ringing artefacts near the edges are present in
the recovered images[30, 31].

• Constrained Least-Squares Filter

Constrained Least-Squares Filter (CLSF) is another image restoration technique from
the family of Least-Squares Filters. In CLSF the smoothness of the restored image is
maximized subject to a constraint on the fidelity of the restored image. CLSF can
expressed in mathematical as fellow;

 H* 
F 'CLSF =  2 
G (2.11)
 | H | + | C | 
2

HereC is the Fourier Transform (FT) of the high pass filter such Laplacian, H* is the
complex conjugate of H and α is the regularization parameter also known is
Lagrange multiplier. The constraint on the least-squares explained here below.

g − Hf '
2

min imize J ( f ') = (2.12)

f' 2

2
subject to Cf '  (2.13)
2

Eq. 2.12 presents the fidelity of the data and Eq. 2.13 is the prior knowledge. Here the
ε is the threshold on the higher frequency (to keep smoothen the restored image). In
Eq. 2.11 if α= 0 then the CLS filter becomes the inverse filter. The threshold on the
higher frequency should be set before the restoration process [27, 32, 33].

34
2.7.3. ITERATIVE BLIND DECONVOLUTION (IBD) TECHNIQUE

The Iterative Blind Deconvolution (IBD) method makes use of the Fast Fourier Transform (FFT)
and the deterministic constraints in the form of non-negativity and finite support constraints. The
iterative algorithms presented in[30, 34-37] and its’ basics procedure is shown in Fig. 2.9. The
image estimate is denoted by f’, the PSF estimate by h’, and the linearly degraded image by g.
Capital letters represent FFT versions of the corresponding signals. Subscript r denotes the
iteration number of the algorithm.

The iterative process can be summarized as follows:

1. First, a non-negative valued initial estimate f0 is input to the process.

2. This is Fourier transformed to give F’r which is then inverted to form an inverse filter and
used to form a new estimate of G, G’r.
3. G’r is then transformed by inverse Fast Fourier Transform (iFFT) to give gr.
4. Then image non-negativity constraints are imposed to reveal a positive constrained
estimate of gr.
5. gr is then Fourier transformed to give the spectrum of G’r.
6. The spectrum G’r is then inverted to form an inverse filter and multiplied by H’r to
estimate F’r.
7. F’r is then inverse Fourier transformed to give fr.
8. Image constraints are applied and an estimate of the image f’r is achieved, this completes
single iteration of the algorithm.

The iterative loop is repeated until two positive functions with the required convolution have
been found.

Unfortunately, the IBD algorithm suffers from two main problems:

• The inverse filter is difficult to define in regions where the inverted function possesses
regions with low values.
• Spectral zeros at frequencies in either F’rorG’r provide no information about that spatial
frequency being a part of the blurring process.

35
 Fig.0.8: Block diagram of Iterative Blind Deconvolution (IBD) algorithm [6, 7]

Implementation of this basic algorithm differs on the assumption on the true image and PSF,
implementation of the assumptions and application in mind[35, 38]. The IBD method is popular
because of its low complexity[34, 39]. Another advantage of this technique is its robustness to
noise which results from the ill-posed nature of the blind image deconvolution problem. IBD
algorithm also suffers from uncertain uniqueness, convergence, instability and sensitivity to
initial image estimate[35].

2.7.4. RICHARDSON LUCY ALGORITHM

In this algorithm, the images (both original and degraded) and PSF (blurring function) are treated
as probability-frequency functions. It is an iterative deblurring technique based on Bays’

36
theorem. In Eq. 2.3 F (original image), H (PSF) and G (restored image), by applying Bays’
theorem as

P(Gk | F ) P( F )
P( F | Gk ) =
i, j P(Gk | F ) P( F ) (2.14)

Eq. 2.14 gives initial estimate of the recovered image, which suggests an iterative restoration
method given below.

P( F | Gk ) P(Gk )
Pk +1 ( F ) = Pk ( F )
k i, j P( F | Gk ) P( F ) (2.15)

Richardson Lucy algorithm requires initial gauss for the blurring kernel[40]. Lesser iterations
result in the presence of residual blur while more iterations result in the corruption of deblurred
image.

2.7.5. REGULARIZATION METHOD

Looking into the convolution model of blurring presented in Eq. 2.3, the image estimate through
inverse filtering is given by Eq. 2.16as follows;

G V
F'= =F+ (2.16)
H H

The restoration error for this model is given by Eq. 2.16,

2
V V
F '− F = = (2.17)
H H

Due to the ill-posed inverse problem, the restoration error will take very large values,
particularly amplifying the high frequency noise [41]. Due to this problem, the system defined in
Eq. 2.16 yields solutions at points where the amplified high frequency noise masks the desired
solution F. In accordance with the regularization theory [42]. If prior information about the noise
or original data can be incorporated, physically meaningful solutions to the ill-posed problem can
be achieved. If we represent a regularization operator in the frequency domain by L then the

37
regularized restored image can be obtained by minimizing the cost function in Eq. 2.17 in
accordance with [33].

LF ' (2.18)

Provided the norm of the noise is known, the term is given in Eq. 2.19

G − HF ' = V (2.19)

The above constrained minimization problem leads to the constrained least-squares or Tikhonov-
Miller regularized solution[33, 42].

2.7.6. TOTAL VARIATION ALGORITHM

This algorithm is based on numerical algorithm for solving total variation constrained least-
squares problems. The concept of the algorithm is based on an augmented Lagrangian method
proposed in[43], and it is a variation of the popularly known Alternating Direction Methods of
Multipliers (ADMM) [43, 44]. In particular, this technique solves the following minimization
problems.

u
Hf − g
2
minimize ( f ) 2
+ f TV (2.20)
2

minimize ( f ) u Hf − g 1 + f TV (2.21)

H is a circulant matrix denoting a spatially invariant linear operator, μ is a regularization

parameter, and || f ||TV is the total variation norm of the data f, defined as;

f TV
=   x [ Dx f ]2k +  y [ Dy f ]2k +  t [ Dt f ]2k (2.22)
k

f k is the k-th entry of the vector f. The operators D x , Dy and Dt are the gradient operators along
the horizontal, vertical and temporal directions. The relative emphasis of Dx , Dy and Dt can be
controlled by  x ,  y and  t respectively.

38
3. OVERVIEW OF OPTIMIZATION ALGORITHMS
Digital Image Restoration is different operations applying on images to get the information or to
achieve the clean approximation of actual image which is lost due to noise, miss-focus and blur
or motion blur. Until now, the PSF for uniform blurring has been estimated using the functional
form for some common types of blur i.e. the Gaussian blur, motion blur and out-of-focus blur. In
real life, especially in the case of motion blur resulting from camera handshake, the blur follows
convoluted paths resulting in complex shaped PSFs.

There exist various schemes in the literature dealing with the restoration of such images with
some listed in [2, 9, 45, 46]. Their shape cannot be easily modelled by a simple equation or
defined by a mathematical model for a set of its parameter(s). PSFs in such cases exhibit an
arbitrary shape and deblurring in this case requires intricate ways to estimate the blur and
effectively remove it. Arbitrarily shaped PSFs have been shown to exist in the case of
atmospheric turbulence blur as well [20].

The focus of the study to investigate into methods which can help in the quick convergence of
the blind image restoration scheme visualized for arbitrarily shaped PSF estimation scheme.

Different type of work has been done on image restoration as B. DEEVENA RAJU et al. [47] in
their work presents that previous image restoration technique restores the cracks of images by
exploiting the discrete wavelet transform (DWT) and Particle swarm optimization (PSO)
methods.

This image restoration technique presented previously is based on PSO-DWT restored the
images with best quality, but the peak signal to noise ratio (PSNR) value of the restored images
and the threshold selection by PSO reduced the PSO-DWT technique performance.To avoid
these problems, Improved Particle Swarm Optimization (IPSO) algorithm technique is proposed
in their work. It’s a new image restoration technique which is proposed with DWT and IPSO
techniques as figure 1 and 2 shows the restored images with IPSO-DWT and PSO-DWT.

In this technique, the given images are restored by the DWT and the thresholds values are
selected by the Improved Particle Swarm Optimization (IPSO) method. Also, figure 3 shows
original Image and the result restored Image with PSO-DWT and IPSO-DWT respectively.

39
The IPSO technique in theirresearch work solved the PSO technique problems by considering
the bad components in the computational process. The technique is implemented and the results
are analyzed completely with their PSNR values. This shows good and effective results of the
image restored technique in restoring the images and achieving improvement in PSNR value as
Figure 3 is the Comparison of PSO-DWT and IPSO-DWT. PSO also used for motion deblurring
of single photograph [19, 46].Also, the performance of the IPSO-DWT technique is found by
comparing with already existing image restoration PSO-DWT technique.

Fig. 3.1: Input image and restored images with IPSO-DWT and PSO-DWT

40
 (a) (b)(c)
Fig.3.2: (a) Original Image (b) Restored Image with PSO-DWT (c) Restored Image with IPSO-
DWT

Fig. 3.3: Comparison of PSO-DWT and IPSO-DWT

41
Nimali Rajakaruna et al. work presents in[48] that a heuristic technique namely Particle Swarm
Optimization is developed to deblurred the obtained image be reversing its motion effect. Most
of the Smartphone’s contains 3 axis accelerometers and gyroscopes. The Particle Swarm
Optimization uses the data of accelerometers and gyroscopes to find out the motion vector of the
image, which is calculated due to the motion of the Smartphone during the time of image
capturing.
Their research work results show that deblurring is successfully done using motion vector which
is derived using the data of accelerometers and gyroscopes. And that deblurred images are used
without any difficulty for the path identification in vision based navigation. Figure 4 shows
different deblurring techniques for a single photograph.

Also, these images used for indoor/outdoor blind and vision navigation[49]. The proposed
method in their work also compared with others deblurring techniques, results shows image
quality was better in proposed method as Figure 5 shows different deblurring techniques and
Improvements.

(a) (b) (c) (d)

(e) (f)

Fig.3.4: Images with different deblurring techniques (a) Original Image (b) Blind deconvolution (c)
Wiener filter (d) Lucy-Richardson (e) Regularized filter (f) Proposed method.

42
 Fig.3.5: Different deblurring techniques and Improvements

In [50] introduced Ant Colony Optimization (ACO) algorithm first time. Their algorithm based
on artificial ants which find out the shortest feasible path by using the amount of pheromone
hormones deposited on the edges of the graph as Figure 3.6 shows the optimal solution of the
problem. Computer results of their shows that the ant colony algorithm is capable of getting
good solutions to both symmetric and asymmetric instances of the problem. The method is an
example, like simulated annealing, neural networks, and evolutionary computation, of the
successful use of a natural metaphor to design an optimization algorithm.

Fig.3.6: Graph shows the CPU time for optimal solution. Continues solid line shows absence of the
algorithm to the problem. Dot line shows the algorithm optimal solution

43
There are Various heuristic optimization schemes used for different types of optimization
problems listed as:

• Simulated annealing (SA)

• Tabu search (TS)
• Differential evolution (DE)
• Memetic algorithms (MA)
• Particle swarm optimization (PSO)
• Graduated optimization (GO)
• Ant Colony Optimization (ACO)
• Artificial Neural Networks (ANNs)
• Support Vector Machines (SVMs)
• Genetic algorithm (GA)
• Bees Algorithm (BA)

This is not the end, apart from these other optimizations schemes also exist. Now the question is
why we select these three optimization schemes from above list. The reason was, after studying
all the above schemes we come to know that these three optimization schemes ACO, GA and
PSO have low coding complexity, easily to understand and faster as compared to other schemes.

44
4. OPTIMIZATION SCHEMES
The three optimization techniques namely GA, PSO and ACO are presented briefly as follows:
4.1. GENETIC ALGORITHM (GA)

Genetic Algorithm is a heuristic search optimization algorithm actually based on the ideas of
genetics and natural selection. The algorithm was first proposed by John Holland, some students
and colleagues for the optimization problems. Since been GA is using for large number of
optimization problems. The design of GA is based on the population which is all the possible
solutions to the required task. Figure 1 shows the GA flowchart. The idea and technique of the
GA is completely based on the natural genetics, in which the pairs combined to produced new
offspring.

In offspring, each individual got the optimal value, the individuals with the most optimal value
further paired to make the results more optimal and make their offspring fitter like a Darwin
theory of evolution “survival to the fittest”. The whole process repeated generation to generation
to make the results more optimal and acceptable [51]. In GA, the population is considered as
chromosomes found in DNA transfer from generation to generation. Individuals in offspring
shows point in a search space and the optimal solution to it.

45
 Initialize Population

Calculate fitness Function

NO Select the fittest chromosomes

Generate Child Population

(Crossover / Mutation)

Stopping Criteria Fulfilled

YES

Terminate with Best Result

Fig. 4.1: GA stages of execution.

• The GA considered individuals with most optimal value for offspring

• In process, some offspring has better value than their parent and thus chose for further
pairing
• By selecting the individuals with best values makes the environment more favourable for the
solution to the problem
• Thus, only those individuals are promoted which are the fittest

4.2. ANT COLONY OPTIMIZATION (ACO)

Marco Dorigo in 1992 firstly presented ACO in his PhD research work in [50]. They introduce
first time the ant colony algorithm which was capable to solve the travelling salesman problem.
The algorithm based on artificial ants which find out the shortest feasible path by using the
amount of pheromone hormones deposited on the edges which shows the optimal solution of the
problem [52].

46
 Begin

Initialize Ants

Search Solution

NO Pheromone update

Addition actions

Exit Condition

YES

End

Fig. 4.2: ACO stages of execution

At the start, each ant was put on the city randomly. To achieve a feasible solution, ants use
probabilistic decision rule to select the city to be visited[52]. When an ant k states in city i and
constructs the partial solution, the probability of ant moving to the next city is given by

(4.1)

τij is the intensity of trails between paths (i.j) and ηij is the visibility of paths (i, j), and ηij=1/dij.
Jk(i) is a set of cities which have to be visited when the ant is at city i. α and β respectively are

47
two positive parameters that control the relative cost of the pheromone trail and of the visibility.
The pheromone amount on each path will be found with equation.

(4.2)

The aim of algorithm was to find out most favourable path in a graph. This was actually based on
the social behaviour of ants finding the path in their colonies and a source of food. Ants are
naturally trained to find out the shortest path to bring food to their nest. They communicate
through a special chemical named pheromone released through their bodies. Pheromone is
released when they leave their nests and when they find food, so that other ants start visiting the
food source by following the pheromone trail [49].

The proposed algorithm for travelling salesman problem is described as:

InitiallySet parameters and initialize pheromone

Tofind out the maximum entropy

Iterativeloops used here to check the whole info repeatedly

Every node has one ant so we have to position it as to distribution strategy

The For loop will go from 1 to m

At the start each ant should go through different path

Repeat

Select next node j to be visited

Until the ant complete the tour

Terminate the For Loop

Update parameters

End

48
 (a)

(b)

(c)

(d)

Fig.4.3. (a) Ants find out the path between nest and food source. (b) when obstacle comes in their
path: then Ants decide by their intelligent behavior whether to turn left or right. (c) More ants pass
through the shortest path so the pheromone amount is stronger than other paths. (d) Ant follow the
path with much amount of pheromone and by this way they choose the shortest path [1]

49
The ants follow the path containing pheromones; they move according to the number of
pheromones depending on the amount of the chemicals they follow the path containing more
pheromones. In turn, the shortest path followed by the majority of ants has a higher level of
pheromone residue and has a higher probability of being selected thus leading the algorithm
towards a quicker convergence. ACO has been used by computer scientists to solve many
numerical problems at that time and adopted broadly at engineering level to emerge the problems
[49].

4.3. PARTICLE SWARM OPTIMIZATION (PSO)

Particle Swarm Optimization (PSO) is a heuristic optimization algorithm, proposed by Kennedy

and Eberhart in 1995 [53]. Now a day it’s the most favourable algorithm for the optimization
techniques. The algorithm was first used for the simulation of social behaviour, for
representation of the movement of bird flock or fish school [54]. Later on, the algorithm was
more simplified and was considered performing optimization. The algorithm is based on swarm
intelligence and it is based on the social behaviour of the bird and fish flock. The birds have their
own social behaviour to find the food at different places, either they scattered completely or go
together in search of food that where they can find. The birds move to search the food from one
place to another, there is always a bird that can better smell the food and they are perceptible of
the different place that contains the food. According to their natural behaviour they always
transmit the information, so at the time of searching the food they always share good information
any time. The concern of this social behaviour to Particle Swarm Optimization algorithm is that
[55]

• The particle swarm is considered as the bird swarm

• Initially Solution development is considered as the bird food searching

• Most favourable information is considered as the best solution

• And the food sources of the bird are considered as the most optimist solution during the
whole process

The particles in PSO cooperate with each other to find out the most optimist solution. Also this
algorithm used for most complex optimization problem [55]. As PSO is an optimization
algorithm and used for most complex optimistic problems, large research has been done to study

50
and improve the performance of the algorithm. All the studies focus to understand better the PSO
and its proper working and how to improve the algorithm. Study shows that PSO has its own
parameters which are inertia weight, velocity, particle swarm size and acceleration coefficients.

Initialize particles

Calculate fitness values

Assign best particles values

Calculate velocity for each

NO particle

Use each particle velocity

value to update its data value

Maximum reached?

YES

End

Fig. 4.4: PSO stages of execution

51
5. EXPERIMENTAL RESULTS FOR OPTIMIZATION OF
PSF ESTIMATION
The aim of this research is to investigate into methods which can aid in the faster convergence of
the blind image restoration scheme. The chapter is organized as follows: Section 5.1 presents to
the reader the experimental results and setups of Ant Colony Optimization (ACO), Genetic
Algorithm (GA) and Particle Swarm Optimization (PSO) which have been utilized in this
research work.

5.1. EXPERIMENTAL SETUP AND RESULTS

In this research work, we compared the proposed algorithms ACO algorithm, PSO and GA
results with each other. The results clearly show that PSO is better in performance than ACO and
Genetic algorithm, it gives better results in minimum time with different iteration. The
experiment shows that PSO in this research work attained better results for single image
deblurring, and the efficiency of solution is higher than ACO and GA, also the convergence
speed is better than that of other algorithms.

Figures 5.1, 5.2, 5.3 and 5.4 shows the blind deblurring restoration for the image using GA, ACO
and PSO estimated PSFs. The results for PSO based estimation appears high quality as compared
to the ones obtained by ACO and GA. The estimates of ACO and GA suffer from ringing around
the border of the images.
(a) (b)

(a) (b) (c) (d)

(c) (d) (e) (f)

52
 Fig. 5.1: (a)Real motion blurred image(b) Deblurred using GA (c) Deblurred using ACO
(d)Deblurred using PSO scheme.

Figure shows the real blurred image as:

(a) (b)

(c) (d)
Fig. 5.2: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO

(a)

53
 (b)

(c)

(d)
Fig. 5.3: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO

(a) (b)

54
 (c) (d)
Fig. 5.4: (a) Real blurred Image (b) Deblurred image with GA (c) Deblurred Image with ACO
(d)Deblurred Image with PSO

Table 5.1 presents the convergence rates for the three optimization schemes. PSO based blind
deblurring scheme is able to converge faster as compared to the GA and PSO schemes. This is
primarily due to the attraction behaviour of PSO such that other members of the swarm are
attracted towards the best solution. PSO is also comparatively easier to implement and has a
lower coding complexity. The research investigates for single image at a time to complete the
whole process.

55
 Blurring Estimated Values
Image
Radius GA ACO PSO
Barbara 12.35 13.25 12.1 12.05
Washsat 9.71 9 10.72 9.15
Boat 15.1 14.73 15.01 15.1
Cameraman 17.45 17.38 16.6 17.67
Goldhill 22.25 21.5 19.85 21.45
Lena 17.3 17.25 18.74 16.54
Mandrill 19.7 19.75 18.9 19.65
Peppers 27.79 27.75 28.4 29.26

Convergence Time (s)

Image
GA ACO PSO
Barbara 43 37 33
Washsat 36 39 32
Boat 49 41 37
Cameraman 23 27 20
Goldhill 32 24 27
Lena 28 22 19
Mandrill 41 33 29
Peppers 30 36 25
Average 35.25 32.37 27.75

Table. 5.1: Results of ACO, GA and PSO

This paper used different heuristic optimization algorithms for single image deblurring to solve
the problem. The main contribution of this paper is a study three heuristic algorithms ACO, PSO
and GA and compare empirically the convergence of these schemes. The experimental result and
performance comparison of the algorithms shows that PSO has much better efficiency than ACO
and genetic algorithm. The implementation of PSO and ACO algorithm is complicated as
compared to Genetic algorithm but the results are efficient. In future, we can compare the results
with more heuristic algorithms for the improvement of performance and betterment. Also, we
will try to implement the whole system in parallel through CUDA and to make it more effective
we will check it on Graphic Processing Unit (GPU).

56
6. DISCUSSION AND ANALYSIS
Deconvolution algorithms cannot apply directly to non-Uniform motion deblurring to find out
the blur kernel. The purpose of the study was to apply different computational algorithm to a
single blurred image to restore and make it informative. Applying different algorithms, we see
the best approach to the problem, we compared the results and observedthat which algorithm is
working fast and efficient. It’s make a good contribution that which optimization technique is
best for the specific problem.
Apart from genetic algorithm the research deals with different optimization scheme for the blind
deblurring of space invariant blur on a single image. A challenging task was to investigate the
blind deblurring of a single image as an input.Different optimization scheme used to approach
the problem to see that which optimization scheme converge quickly for the blind image
restoration scheme visualized for arbitrarily shaped PSF estimation scheme.
PSF estimation techniques are usually slow as the search space is large and they are further
marred by the presence of large image data especially multi-channel images. To reduce the time
taken to quickly approximate the true PSF, parallelization of the algorithm can be done but this is
not encouraging. Rather, optimization techniques are sought which aim at converging quickly
and estimating the PSF fast.

Investigation of the evolutionary algorithms for search optimization shows that PSO converges
fast as compare to GA and ACO.

57
7. CONCLUSIONS AND FUTURE WORK
In this research study, we explored the problem of blind image restoration also known as Blind
Image Deconvolution. The research work focused on the analysis of the existing BID schemes
and to find a way for the optimization of these restoration methods. A brief discussion of the
important findings and the suggestion for the future work are suggested here in this chapter.

7.1. CONCLUSIONS

This research work used different heuristic optimization algorithms for single image deblurring
to solve the problem. The main contribution of this research is the study three heuristic
algorithms ACO, PSO and GA and compare empirically the convergence of these schemes. The
experimental result and performance comparison of the algorithms shows that PSO has much
better efficiency than ACO and genetic algorithm. The implementation of PSO and ACO
algorithm is complicated as compared to Genetic algorithm but the results are efficient. In future,
we can compare the results with more heuristic algorithms for the improvement of performance
and betterment. Also, we will try to implement the whole system in parallel through CUDA and
to make it more effective we will check it on Graphic Processing Unit (GPU).

7.2. RESEARCH FINDINGS

The research contributed to the field of Digital Image Restoration for BID by introducing the
Quick Convergence of PSF Estimation for Single Image Blind Deblurring.
In this work,apart from genetic algorithm various optimization schemes were used for the blind
deblurring of a single image. A challenging task is to investigate the blind deblurring of a single
image as an input. Different optimization scheme was used to approach the problem to see that
which optimization scheme converge quickly for the blind image restoration scheme visualized
for arbitrarily shaped PSF estimation scheme.
The main findings are using various optimization schemes like ACO, GA and PSO for PSF
estimation of the single BID. By comparing all the results of the schemes, we found clearly that
PSO converge the problem faster than that of ACO and GA. Also, the estimated valves of PSO is
much better in most cases as compared to other schemes.

58
7.3. FUTURE WORK

In future, we will extend the work to investigate more optimization schemes for single image
blind blurring. for the purpose to compare the results with more heuristic algorithms for the
improvement of performance and betterment.
This work currently deal with three optimization schemes mainly as ACO, GA and PSO which
was used for the fast convergence for single BID. Results achieved from the schemes were
compared in between and after good observation PSO converge faster as compared to other
schemes. But still three optimization schemes results were compared with each other, we will
pursue the work in near future for further optimization schemes.
Currently this research work deal only with non-parametric PSF. The research does not deal with
the parametric PSF and denoising (noise removal) of the image may it be artificial or deblurring
noise. Also, ringing effects of an image is out of scope for the research. In onward work the
circle of the research will also deal with the parametric PSF.
We will further try to extend the work to utilize the noise removal or denoising factors. The BID
schemes will also be used then for noisy images. Also, we will try to implement the whole
system in parallel through CUDA and to make it more effective we will check it on Graphic
Processing Unit (GPU).

59
REFERENCES

1. Kang, W., et al., Fast digital zooming system using directionally adaptive image interpolation
and restoration. Springerplus, 2014. 3: p. 9.

2. Whyte, O., et al., Non-uniform Deblurring for Shaken Images. International Journal of Computer
Vision, 2012. 98(2): p. 168-186.

3. Shan, Q., J. Jia, and A. Agarwala, High-Quality Motion Deblurring From a Single Image. ACM
Transactions on Graphics, 2008. 27(3): p. 10.

4. Cruz-Barbosa, R. and A. Vellido, Semi-Supervised Analysis of Human Brain Tumours from

Partially Labelled MRS Information, Using Manifold Learning Models. International Journal of
Neural Systems, 2011. 21(1): p. 17-29.

5. Glover, G.H., Deconvolution of Impulse Response in Event-related BOLD fMRI. Neuroimage,

1999. 9(4): p. 416-429.

6. Taxt, T. and J. Strand, Two-Dimensional Noise-Robust Blind Deconvolution of Ultrasound

Images. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control, 2001. 48(4): p.
861-866.

7. Shi, L.I., Z. Bao, and S.U.N. Hui, Restoration of Aerial Multiple Blurred Images. Optics and
Precision Engineering, 2009. 17(5): p. 1161-1170.

8. Jun, W. and C. Danqing, Deblurring Texture Extraction from Digital Aerial Image by Reforming
a Steep Edge Curve. Geo-spatial Information Science, 2005. 8(1): p. 39-44.

9. Fergus, R., et al., Removing Camera Shake from a Single Photograph. ACM Transactions on
Graphics, 2006. 25(3): p. 787-794.

10. Cho, S. and S. Lee, Fast Motion Deblurring. ACM Transactions on Graphics, 2009. 28(5).

11. Agrawal, A., Y. Xu, and R. Raskar, Invertible Motion Blur in Video. ACM Transactions on
Graphics, 2009. 28(3).

12. Caron, J.N., Blind Deconvolution of Audio-Frequency Signals Using the Self-Deconvolving Data
Restoration Algorithm (SeDDaRA). Journal of the Acoustical Society of America, 2004. 116(1):
p. 373-378.

13. Ligorria, J.P. and C.J. Ammon, Iterative Deconvolution and Receiver-Function Estimation.
Bulletin of the Seismological Society of America, 1999. 89(5): p. 1395-1400.

60
14. Ramsey, M.S. and P.R. Christensen, Mineral Abundance Determination: Quantitative
Deconvolution of Thermal Emission Spectra. Journal of Geophysical Research-Solid Earth, 1998.
103(B1): p. 577-596.

15. Jefferies, S.M., et al., Blind Deconvolution in Optical Diffusion Tomography. Optics Express,
2002. 10(1): p. 46-53.

16. Jefferies, S.M. and J.C. Christou, Restoration of Astronomical Images by Iterative Blind
Deconvolution. Astrophysical Journal, 1993. 415(2): p. 862-874.

17. Desidera, G., et al., Application of Iterative Blind Deconvolution to the Reconstruction of LBT
LINC-NIRVANA Images. Astronomy & Astrophysics, 2006. 452(2): p. 727-734.

18. 12 August 2017; Available from:

https://en.wikipedia.org/wiki/List_of_government_space_agencies.

19. Song, C.H., et al., PSO Based Motion Deblurring for Single Image. Gecco-2011: Proceedings of
the 13th Annual Genetic and Evolutionary Computation Conference, 2011: p. 85-92.

20. Khan, A., Efficient Methodologies for Single-Image Blind Deconvolution and Deblurring, in
School of Electrical and Electronic Engineering. 2014, The University of Manchester.

21. Somasundaram, K. and P.A. Kalaividya. Brain portion segmentation from Magnetic Resonance
Images(MRI) of human head scan using Richardson Lucy deconvolution and intensity
thresholding. in 2016 International Computer Science and Engineering Conference (ICSEC).
2016.

22. Yang, H., et al., Image deblurring using empirical Wiener filter in the curvelet domain and joint
non-local means filter in the spatial domain. The Imaging Science Journal, 2014. 62(3): p. 178-
185.

23. Senshiki, H., et al. PSF estimation using total variation regularization and shock filter for blind
deconvolution. in 2017 IEEE International Conference on Consumer Electronics (ICCE). 2017.

24. Lee, J.Y., K.W. Ko, and S. Lee, Implementation of an Image Restoration with Block Iteration
Method for Spatially Variant Blur Models, in Advances in Computer Science and Ubiquitous
Computing: CSA-CUTE2016, J.J. Park, et al., Editors. 2017, Springer Singapore: Singapore. p.
414-420.

25. Banham, M.R. and A.K. Katsaggelos, Digital Image Restoration. IEEE Signal Processing
Magazine, 1997. 14(2): p. 24-41.

61
26. Li, D.L. and S. Simske, Atmospheric Turbulence Degraded-Image Restoration by Kurtosis
Minimization. IEEE Geoscience and Remote Sensing Letters, 2009. 6(2): p. 244-247.

27. Gonzalez, R.C. and R.E. Woods, Digital Image Processing. Prentice-Hall Inc, 2002.

28. Lagendijk, R.L. and J. Biemond, Basic Methods for Image Restoration and Identification, in The
Essential Guide to Image Processing. 2009, Academic Press USA: San Diego. p. 323-348.

29. Hufnagel, R.E. and N.R. Stanley, Modulation Transfer Function Associated With Image
Transmission through Turbulent Media. Journal of the Optical Society of America, 1964. 54(1):
p. 52-&.

30. Kundur, D. and D. Hatzinakos, Blind Image Deconvolution. IEEE Signal Processing Magazine,
1996. 13(3): p. 43-64.

31. Stockham Jr, T.G., T.M. Cannon, and R.B. Ingebretsen, Blind deconvolution through digital
signal processing. Proceedings of the IEEE, 1975. 63(4): p. 678-692.

32. Bilgen, M. and H.-S. Hung, Constrained least-squares filtering for noisy images blurred by
random point spread function. Optical Engineering, 1994. 33(6): p. 2020-2023.

33. Hunt, B.R., The Application of Constrained Least-Squares Estimation to Image Restoration by
Digital Computer. IEEE Transactions on Computer, 1973. C-22: p. 805-812.

34. Kundur, D. and D. Hatzinakos, A Novel Blind Deconvolution Scheme for Image Restoration
Using Recursive Filtering. IEEE Transactions on Signal Processing, 1998. 46(2): p. 375-390.

35. Ayers, G.R. and J.C. Dainty, Iterative Blind Deconvolution Method and Its Applications. Optics
Letter, 1988. 13(7): p. 43-64.

36. Biemond, J., R.L. Lagendijk, and R.M. Mersereau, Iterative Methods for Image Deblurring.
Proceedings of the IEEE, 1990. 78(5): p. 856-883.

37. Katsaggelos, A.K., Iterative Image Restoration Algorithms. Optical Engineering, 1989. 28(7): p.
735-748.

38. Davey, B.L.K., R.G. Lane, and R.H.T. Bates, Blind Deconvolution of Noisy Complex-Valued
Image. Optics Communications, 1989. 69(5-6): p. 353-356.

39. Kundur, D. and D. Hatzinakos. A novel recursive filtering method for blind image restoration. in
Proc. IASTED Int. Conf. on Signal and Image Processing. 1995.

40. Richardson, W.H., Bayesian-Based Iterative Method of Image Restoration. Journal of the Optical
Society of America, 1972. 62(1): p. 55-59.

62
41. Lagendijk, R.L., J. Biemond, and D.E. Boekee, Regularized Iterative Image-Restoration with
Ringing Reduction. IEEE Transactions on Acoustics Speech and Signal Processing, 1988. 36(12):
p. 1874-1888.

42. Tikhonov, A.N. and V.Y. Arsenin, Solutions of Ill-Posed Problems. 1977, New York: Wiley.

43. Chan, S.H., et al., An augmented Lagrangian method for total variation video restoration. Image
Processing, IEEE Transactions on, 2011. 20(11): p. 3097-3111.

44. Han, D., X. Yuan, and W. Zhang, An augmented Lagrangian based parallel splitting method for
separable convex minimization with applications to image processing. Mathematics of
Computation, 2014. 83(289): p. 2263-2291.

45. Kim, J.B. and H.J. Kim, GA-based image restoration by isophote constraint optimization. Eurasip
Journal on Applied Signal Processing, 2003. 2003(3): p. 238-243.

46. Jing, W., et al, Motion deblurring for single photograph based on particle swarm optimization.
International Journal of Computational Intelligence Systems: p. p. 1-11.

47. RAJU, B.D.P., P.; PRASAD, G. S, A new Image Reconstruction Technique with Aid of IPSO
(Improved Particle Swarm Optimization) and DWT (Discrete Wavelet Transform). Journal of
Theoretical & Applied Information Technology.

48. Rajakaruna, N., et al., Image Deblurring for Navigation Systems of Vision Impaired People Using
Sensor Fusion Data. 2014 Ieee Ninth International Conference on Intelligent Sensors, Sensor
Networks and Information Processing (Ieee Issnip 2014), 2014.

49. Zhang, Y.D., S.H. Wang, and G.L. Ji, A Comprehensive Survey on Particle Swarm Optimization
Algorithm and Its Applications. Mathematical Problems in Engineering, 2015: p. 38.

50. Dorigo, M. and C. Blum, Ant colony optimization theory: A survey. Theoretical Computer
Science, 2005. 344(2-3): p. 243-278.

51. Man, K.F., K.S. Tang, and S. Kwong, Genetic algorithms: Concepts and applications. Ieee
Transactions on Industrial Electronics, 1996. 43(5): p. 519-534.

52. Wei, X., L. Han, and L. Hong, A Modified Ant Colony Algorithm for Traveling Salesman
Problem. International Journal of Computers Communications & Control, 2014. 9(5): p. 633-643.

53. X.H. Shi, Y.C.L., H.P.Lee b,,C.Lu ,Q.X.Wang Particle swarm optimization-based algorithms for
TSP and generalized TSP. information processing letters, 2007: p. 169-176.

63
54. Bai, Q. and Analysis of Particle Swarm Optimization Algorithm. computer and information
science, 2010. 3(1): p. 180-184.

55. Kao, C.-C., Application of Particle Swarm Optimization in Mechanical Design. 2009: p. 93-116.

64
APPENDICES

A. Appendix A

In this research three optimization schemes are used.

Optimization schemes Code

GA

function varargout = tsp_ga(varargin)

clc
clear all
close all

tstart=tic;
cput=cputime;
load xy.mat
defaultConfig.xy =xy;

defaultConfig.xy = 10*rand(50,2);
defaultConfig.dmat = [];
defaultConfig.popSize = 100;
defaultConfig.numIter = 1e4;
defaultConfig.showProg = true;
defaultConfig.showResult = true;
defaultConfig.showWaitbar = false;

if ~nargin
userConfig = struct();
elseif isstruct(varargin{1})
userConfig = varargin{1};
else
try
userConfig = struct(varargin{:});
catch
error('Expected inputs are either a structure or parameter/value
pairs');
end
end

configStruct = get_config(defaultConfig,userConfig);
xy = configStruct.xy;
dmat = configStruct.dmat;
popSize = configStruct.popSize;
numIter = configStruct.numIter;
showProg = configStruct.showProg;
showResult = configStruct.showResult;
showWaitbar = configStruct.showWaitbar;
if isempty(dmat)
nPoints = size(xy,1);
a = meshgrid(1:nPoints);
dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),nPoints,nPoints);
end

65
 [N,dims] = size(xy);
[nr,nc] = size(dmat);
if N ~= nr || N ~= nc
error('Invalid XY or DMAT inputs!')
end
n = N;
popSize = 4*ceil(popSize/4);
numIter = max(1,round(real(numIter(1))));
showProg = logical(showProg(1));
showResult = logical(showResult(1));
showWaitbar = logical(showWaitbar(1));

pop = zeros(popSize,n);
pop(1,:) = (1:n);
for k = 2:popSize
pop(k,:) = randperm(n);
end

globalMin = Inf;
totalDist = zeros(1,popSize);
distHistory = zeros(1,numIter);
tmpPop = zeros(4,n);
newPop = zeros(popSize,n);
if showProg
figure('Name','TSP_GA | Current Best Solution','Numbertitle','off');
hAx = gca;
end
if showWaitbar
hWait = waitbar(0,'Searching for near-optimal solution ...');
end

for iter = 1:numIter

for p = 1:popSize
d = dmat(pop(p,n),pop(p,1));
for k = 2:n
d = d + dmat(pop(p,k-1),pop(p,k));
end
totalDist(p) = d;
end

[minDist,index] = min(totalDist);
distHistory(iter) = minDist;
if minDist < globalMin
globalMin = minDist;
optRoute = pop(index,:);
if showProg
rte = optRoute([1:n 1]);
if dims > 2, plot3(hAx,xy(rte,1),xy(rte,2),xy(rte,3),'r.-');
else plot(hAx,xy(rte,1),xy(rte,2),'r.-'); end
title(hAx,sprintf('Total Distance = %1.4f, Iteration =
%d',minDist,iter));
drawnow;
end
end

66
 randomOrder = randperm(popSize);
for p = 4:4:popSize
rtes = pop(randomOrder(p-3:p),:);
dists = totalDist(randomOrder(p-3:p));
[ignore,idx] = min(dists);
bestOf4Route = rtes(idx,:);
routeInsertionPoints = sort(ceil(n*rand(1,2)));
I = routeInsertionPoints(1);
J = routeInsertionPoints(2);
for k = 1:4
tmpPop(k,:) = bestOf4Route;
switch k
case 2
tmpPop(k,I:J) = tmpPop(k,J:-1:I);
case 3
tmpPop(k,[I J]) = tmpPop(k,[J I]);
case 4
tmpPop(k,I:J) = tmpPop(k,[I+1:J I]);
otherwise
end
end
newPop(p-3:p,:) = tmpPop;
end
pop = newPop;

if showWaitbar && ~mod(iter,ceil(numIter/325))

waitbar(iter/numIter,hWait);
end

end

telapsed=toc(tstart);
[telapsed cputime-cput]

if showWaitbar
close(hWait);
end
if showResult
figure('Name','TSP_GA | Results','Numbertitle','off');
subplot(2,2,1);
pclr = ~get(0,'DefaultAxesColor');
if dims > 2, plot3(xy(:,1),xy(:,2),xy(:,3),'.','Color',pclr);
else plot(xy(:,1),xy(:,2),'.','Color',pclr);
end

title('City Locations');
subplot(2,2,2);
imagesc(dmat(optRoute,optRoute));
title('Distance Matrix');
subplot(2,2,3);
rte = optRoute([1:n 1]);
if dims > 2, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'r.-');

67
else plot(xy(rte,1),xy(rte,2),'r.-'); end
title(sprintf('Total Distance = %1.4f',minDist));
subplot(2,2,4);
plot(distHistory,'b','LineWidth',2);
title('Best Solution History');
set(gca,'XLim',[0 numIter+1],'YLim',[0 1.1*max([1 distHistory])]);
end

% Return Output
if nargout
resultStruct = struct( ...
'xy', xy, ...
'dmat', dmat, ...
'popSize', popSize, ...
'numIter', numIter, ...
'showProg', showProg, ...
'showResult', showResult, ...
'showWaitbar', showWaitbar, ...
'optRoute', optRoute, ...
'minDist', minDist);

varargout = {resultStruct};
end

end

function config = get_config(defaultConfig,userConfig)

config = defaultConfig;

defaultFields = fieldnames(defaultConfig);

userFields = fieldnames(userConfig);
nUserFields = length(userFields);

for i = 1:nUserFields
userField = userFields{i};
isField = strcmpi(defaultFields,userField);
if nnz(isField) == 1
thisField = defaultFields{isField};
config.(thisField) = userConfig.(userField);
end
end

end

ACO

clc
clear all
close all

tstart=tic;
cput=cputime;

68
 [x,y,d,t,h,iter,alpha,beta,e,m,n,el]=ants_information;
for i=1:iter
[app]=ants_primaryplacing(m,n);
[at]=ants_cycle(app,m,n,h,t,alpha,beta);
at=horzcat(at,at(:,1));
[cost,f]=ants_cost(m,n,d,at,el);
[t]=ants_traceupdating(m,n,t,at,f,e);
costoa(i)=mean(cost);
[mincost(i),number]=min(cost);besttour(i,:)=at(number,:);
iteration(i)=i;
end

telapsed=toc(tstart);
[telapsed cputime-cput]

subplot(2,1,1);plot(iteration,costoa);
title('average of cost (distance) versus number of cycles');
xlabel('iteration');
ylabel('distance');
[k,l]=min(mincost);
for i=1:n+1
X(i)=x(besttour(l,i));
Y(i)=y(besttour(l,i));
end
subplot(2,1,2);plot(X,Y,'--rs','LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',10)
xlabel('X');ylabel('y');axis('equal');

figure(2)
plot(X,Y,'--rs')

xlabel('X');ylabel('y');axis('equal');

for i=1:n
text(X(i)+.5,Y(i),['\leftarrow node ',num2str(besttour(l,i))]);
end
title(['optimum course by the length of ',num2str(k)]);

function [cost,f]=ants_cost(m,n,d,at,el);
for i=1:m
s=0;
for j=1:n
s=s+d(at(i,j),at(i,j+1));
end
f(i)=s;
end
cost=f;
f=f-el*min(f);
function [at]=ants_cycle(app,m,n,h,t,alpha,beta);
for i=1:m
mh=h;
for j=1:n-1

69
 c=app(i,j);
mh(:,c)=0;
temp=(t(c,:).^beta).*(mh(c,:).^alpha);
s=(sum(temp));
p=(1/s).*temp;
r=rand;
s=0;
for k=1:n
s=s+p(k);
if r<=s
app(i,j+1)=k;
break
end
end
end
end
at=app;

function [x,y,d,t,h,iter,alpha,beta,e,m,n,el]=ants_information;

iter=100;
m=100;

load xy.mat
x=xy(:,1)';
y=xy(:,2)';

n=length(x);
for i=1:n
for j=1:n
d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
end
end
e=.1;
alpha=1;
beta=5;
for i=1:n
for j=1:n
if d(i,j)==0
h(i,j)=0;
else
h(i,j)=1/d(i,j);
end
end
end
t=0.0001*ones(n);
el=.96;
End

function [app]=ants_primaryplacing(m,n);
rand('state',sum(100*clock));
for i=1:m
app(i,1)=fix(1+rand*(n-1));
end

70
function [t]=ants_traceupdating(m,n,t,at,f,e);
for i=1:m
for j=1:n
dt=1/f(i);
t(at(i,j),at(i,j+1))=(1-e)*t(at(i,j),at(i,j+1))+dt;
end
end

PSO

function A = ackley(x)
c1 = 20;
c2 = 0.2;
c3 = 2*pi;
A = -c1 * exp(-c2*sqrt(mean(x.^2))) - exp(mean(cos(c3*x))) + c1 + exp(1);
end

function path = applyEdgeExchange( exchList, path )

if isempty(path) || isempty(exchList)
return
end

l = 1;
if size(exchList,1) > 1
while true

if exchList(l,:) == exchList(l+1,:)
exchList(l+1,:) = [];
end

if l+1 >= size(exchList,1)

break
else
l = l + 1;
end
end
end

for e=1:1:size(exchList,1)
path = edgeR( path, exchList(e,1), exchList(e,2) );
end
end

matlabpool('local', 4);

load('data/datasets.mat');

test = [0 0 0 0];

71
parfor l=1:1:4
for i=1:1:1000
test(1,l) = test(1,l) + 1;
end
end

save('results/test.mat','test');

clearvars

matlabpool close

function result = change_in_path_length(path, i, j, distances)

before = distances(path(r(i, path)), path(r(i+1,path))) +

distances(path(r(i+1+j, path)), path(r(i+2+j, path)));

after = distances(path(r(i, path)), path(r(i+1+j, path))) +

distances(path(r(i+1, path)), path(r(i+2+j, path)));

result = after - before;

end

function path = change_path(path, i, j)

oldPath = path;

path(r(i+1, path)) = oldPath(r(i+1+j, path));

path(r(i+1+j, path)) = oldPath(r(i+1, path));

for k=1:j-1
path(r(i+1+k, path)) = oldPath(r(i+1+j-k, path));
end
end

function d = distance(a)
n = length(a);
d = zeros(n,n);
for i = 1:n
for j = i+1:n
d(i,j)=pdist([a(i,:);a(j,:)],'euclidean');
d(j,i)=d(i,j);
end;
end
end

function dist = distancePath(distances, path)

dist = 0.0;

for e=1:1:size(path,2)-1

dist = dist + distances(path(e), path(e+1));

end

72
 dist = dist + distances(path(end), path(1));
end

function newPath = edgeR( path, i, j )

if i < j

if i > 1
front = path(1,1:i-1);
else
front = [];
end

list = path(1,i:j);

if j < size(path,2)
back = path(1,j+1:end);
else
back = [];
end

newPath = [front fliplr(list) back];

else
newPath = path;
end

end

e1 = 1;
e2 = 2;

particlePos = [51.5 56];

newVelocity = [5 -3];

data = dataset_tspn2DE16_2;

plotEllipse( data(e1,1), data(e1,2), data(e1,3), data(e1,4) );

plotEllipse( data(e2,1), data(e2,2), data(e2,3), data(e2,4) );

plot(particlePos(1,1),particlePos(1,2), '*'); hold on

plot((particlePos(1,1) + newVelocity(1,1)),(particlePos(1,2) +
newVelocity(1,2)), '*'); hold on

plot([particlePos(1,1) (particlePos(1,1) + newVelocity(1,1))],

[particlePos(1,2), (particlePos(1,2) + newVelocity(1,2))], 'b-
','LineWidth',3); hold on

73
XYproj = getNextBoarderPoint( data(e1,:), particlePos, (particlePos +
newVelocity) );

deltaX = XYproj(1,1) - particlePos(1,1);

deltaY = XYproj(1,2) - particlePos(1,2);

if (deltaX / deltaX) > 0

signX = -1;
else
signX = 1;
end

if (deltaY / deltaY) > 0

signY = -1;
else
signY = 1;
end

particlePosNew(1,1) = XYproj(1,1);
particlePosNew(1,2) = XYproj(1,2);

plot(particlePosNew(1,1),particlePosNew(1,2), 'x','MarkerSize',10);

hold on

load('data/datasets.mat');

fullloops = 2;

resultDist = zeros(fullloops, 2, size(datasetname,2));

resultTime = zeros(fullloops, 3, size(datasetname,2));

resultOverview = zeros(size(datasetname,2), 5);

moveOptionsDPSO = struct('bLoc', 0.2, 'bGlob', 0.2, 'randArt', 'randomTemp',

'vRandType', '2opt', 'vRandIter', 0, 'noChangeIterStop', 100);

moveOptionsPSO = struct('initVelocity', 0.0, 'initInertiaWeight', 0.3,

'noChangeCountTh', 0, 'boundaryhandlingPercentage', 0.0, 'noChangeIterStop',
100);

for f = 1

fprintf('dataset: %s\n',datasetname{f});

data = eval(datasetname{f});

optimalPath = zeros(size(data,1), 2, fullloops);

optimalPathLength = Inf(fullloops,1);

for l=1:1:fullloops

tDpso = tic;

74
 [ path, total_length_dpso ] = psoOptDisc( data , 50, 1000,
moveOptionsDPSO);

dpsoResultTime = toc(tDpso) * 1000;

dpsoResultDist = total_length_dpso;

tPso = tic;
[ path, total_length_pso, travelPoints ] = psoOpt( data, path, 50,
1000, moveOptionsPSO);

psoResultTime = toc(tPso) * 1000;

psoResultDist = total_length_pso;

resultDist(l,:,f) = [total_length_dpso total_length_pso];

resultTime(l,:,f) = [dpsoResultTime psoResultTime (dpsoResultTime +

psoResultTime)];

optimalPath(:,:,l) = travelPoints;
optimalPathLength(l,1) = total_length_pso;
end

resultOverview(f,3) = mean(resultDist(:,2,f));

resultOverview(f,1) = (resultOverview(f,3) - optimalValue(f,1)) /

optimalValue(f,1);

resultOverview(f,2) = max(resultDist(:,2,f));

resultOverview(f,5) = min(resultDist(:,2,f));

resultOverview(f,4) = std(resultDist(:,2,f));

fprintf('relative error: %f\n', resultOverview(f,1));

fprintf('max. value: %f\n', resultOverview(f,2));
fprintf('mean value: %f\n', resultOverview(f,3));
fprintf('standard deviation: %f\n', resultOverview(f,4));
fprintf('best solution: %f\n', resultOverview(f,5));

save(['results/optimalPath_' num2str(f) '.mat'], 'optimalPath',

'optimalPathLength');
end

save('results/resultOverview.mat','resultOverview');

clearvars
variablesshowplot_dpsoshowplot_psofltDpsotPsofullloopstotal_length_dpsototal_
length_psomoveOptionsDPSO

load('data/datasets.mat');

75
fullloops = 4;

resultOverview = zeros(size(datasetname,2), 5);

moveOptionsDPSO = struct('bLoc', 0.2, 'bGlob', 0.2, 'randArt', 'randomTemp',

'vRandType', '2opt', 'vRandIter', 0, 'noChangeIterStop', 200, 'particleIter',
100);

moveOptionsPSO = struct('initVelocity', 0.0, 'omega', 0.3, 'noChangeCountTh',

0, 'boundaryhandlingPercentage', 0.0, 'noChangeIterStop', 1001, 'c1', 1.0,
'c2', 1.0, 'particleIter', 100);

for f = 1

fprintf('dataset: %s\n',datasetname{f});

data = eval(datasetname{f});

[ optimalPath, optimalPathLength, optimalPathOrder, resultDist, resultTime ]

= optimization( fullloops, data, moveOptionsDPSO, moveOptionsPSO );

resultOverview(f,3) = mean(resultDist(:,2));

resultOverview(f,1) = (resultOverview(f,3) - optimalValue(f,1)) /

optimalValue(f,1);

resultOverview(f,2) = max(resultDist(:,2));

resultOverview(f,5) = min(resultDist(:,2));
resultOverview(f,4) = std(resultDist(:,2));

end

clearvars

function [ globalBestPath, globalBestDist ] = findGlobalBestDpso( distances,

particlePos, globalBest )

if nargin == 2
globalBest = particlePos(1,:);
end

globalBestPath = globalBest;
globalBestDist = distancePath(distances, globalBestPath);

for d=1:1:size(particlePos,1)

particlePosDist = distancePath(distances, particlePos(d,:));

if particlePosDist < globalBestDist

globalBestPath = particlePos(d,:);

76
 globalBestDist = particlePosDist;
end
end
end

function [ globalBestCoord, globalBestDist ] = findGlobalBestPso( path,

particlePos, globalBest )

if nargin == 2
globalBest = particlePos(:, :, 1);
end

globalBestCoord = globalBest;
globalBestDist = distancePath(distance(globalBest), path);

for d=1:1:size(particlePos,3)

particlePosDist = distancePath(distance(particlePos(:,:,d)), path);

if particlePosDist < globalBestDist

globalBestCoord = particlePos(:,:,d);
globalBestDist = particlePosDist;
end
end
end

function list = getListOfEdgeExchanges( ispath, shallpath )

list = [];

if size(ispath,2) ~= size(shallpath,2)
return
end

if sum( ispath == shallpath ) == size(shallpath,2)

return
end

while sum( ispath == shallpath ) ~= size(shallpath,2)

for i=1:1:size(shallpath,2)

if ispath(1,i) ~= shallpath(1,i)
excEl = find(ispath==shallpath(1,i));
ispath = edgeR(ispath, i, excEl);
list = [list; i excEl];
end
end
end
end

function borderPoint = getNextBoarderPoint( data, oldPoint, newPoint )

77
 borderPoint = [];

rrx = data(1,3)^2;
rry = data(1,4)^2;

x21 = newPoint(1,1) - oldPoint(1,1);

y21 = newPoint(1,2) - oldPoint(1,2);

x10 = oldPoint(1,1) - data(1,1);

y10 = oldPoint(1,2) - data(1,2);

a = (x21 * x21 / rrx) + (y21 * y21 / rry);

b = (x21 * x10 / rrx) + (y21 * y10 / rry);
c = (x10 * x10 / rrx) + (y10 * y10 / rry);
d = b * b - a * (c-1);

if d >= 0
e = sqrt(d);
u1 = (-b-e) / a;
u2 = (-b+e) / a;

if 0 <= u1 && u1 <= 1

borderPoint = [borderPoint; oldPoint(1,1)+x21*u1
oldPoint(1,2)+y21*u1];
end

if 0 <= u2 && u2 <= 1

borderPoint = [borderPoint; oldPoint(1,1)+x21*u2
oldPoint(1,2)+y21*u2];
end
end
end

function particlePos = initSwarmMemberDpso( data, quantity, type )

switch(type)
case'randomNoDup'
particlePos = randomInitializationNoDuplicate( data, quantity );

case'randomWithDup'
particlePos = randomInitializationWithDuplicate( data, quantity );

case'randomFifNoDup'
particlePos = randomInitializationFifNoDuplicate( data, quantity );

otherwise
particlePos = randomInitializationNoDuplicate( data, quantity );
end
end

function particlePos = randomInitializationNoDuplicate( data, quantity )

particlePos = zeros(quantity, size(data,1));

78
if quantity > factorial(size(data,1))
quantity = factorial(size(data,1));
end

for p=1:1:quantity

while true
newPermutation = randsample(1:size(data,1), size(data,1));

if ~ismember(newPermutation, particlePos, 'rows')

particlePos(p,:) = newPermutation;
break
end
end
end
end

function particlePos = randomInitializationWithDuplicate( data, quantity )

particlePos = zeros(quantity, size(data,1));

for p=1:1:quantity

particlePos(p,:) = randsample(1:size(data,1), size(data,1));

end
end

function particlePos = randomInitializationFifNoDuplicate( data, quantity )

particlePos = zeros(quantity, size(data,1));

if quantity > factorial(size(data,1))

quantity = factorial(size(data,1));
end

for p=1:1:quantity

while true
newPermutation = [1 randsample(2:size(data,1), size(data,1)-1)];

if ~ismember(newPermutation, particlePos, 'rows')

particlePos(p,:) = newPermutation;
break
end
end
end
end

function particlePos = initSwarmMemberPso( data, quantity, type )

switch(type)

79
case'random'
particlePos = randomInitialization( data, quantity );

case'fourparts'
particlePos = fourpartsInitialization( data, quantity );

case'center'
particlePos = centerInitialization( data, quantity );

otherwise
particlePos = randomInitialization( data, quantity );
end
end

function particlePos = randomInitialization( data, quantity )

particlePos = zeros(size(data,1), 2, quantity);

for p=1:1:quantity

for c=1:1:size(particlePos,1)

while true

particlePos(c,1,p) = ((data(c,1) + data(c,3)) - (data(c,1) -

data(c,3))) * rand(1) + (data(c,1) - data(c,3));

particlePos(c,2,p) = ((data(c,2) + data(c,4)) - (data(c,2) -

data(c,4))) * rand(1) + (data(c,2) - data(c,4));

if isPointInEllipse( data(c,:), particlePos(c,:,p) ) == true

break
end
end
end
end
end

function particlePos = fourpartsInitialization( data, quantity )

particlePos = zeros(size(data,1), 2, quantity);

for c=1:1:size(particlePos,1)

if quantity == 1
particlePos = randomInitialization( data, quantity );
elseif quantity > 1

aktParticle = 1;

while aktParticle <= quantity

for s=1:1:4

80
if aktParticle > quantity
break
end

upperBound = 0.0;
lowerBound = 0.0;
leftBound = 0.0;
rightBound = 0.0;

switch(s)
case 1
upperBound = data(c,2) + data(c,4);
lowerBound = data(c,2);
leftBound = data(c,1) - data(c,3);
rightBound = data(c,1);

case 2
upperBound = data(c,2) + data(c,4);
lowerBound = data(c,2);
leftBound = data(c,1);
rightBound = data(c,1) + data(c,3);

case 3
upperBound = data(c,2);
lowerBound = data(c,2) - data(c,4);
leftBound = data(c,1);
rightBound = data(c,1) + data(c,3);

case 4
upperBound = data(c,2);
lowerBound = data(c,2) - data(c,4);
leftBound = data(c,1) - data(c,3);
rightBound = data(c,1);
end

while true
particlePos(c,1,aktParticle) = (rightBound -
leftBound) * rand(1) + leftBound;
particlePos(c,2,aktParticle) = (upperBound -
lowerBound) * rand(1) + lowerBound;

if isPointInEllipse( data(c,:), particlePos(c,:,aktParticle) ) == true

break
end
end

aktParticle = aktParticle + 1;
end
end
end
end
end

function particlePos = centerInitialization( data, quantity )

81
 particlePos = zeros(size(data,1), 2, quantity);

for p=1:1:quantity

for c=1:1:size(particlePos,1)

particlePos(c,:,p) = data(c,1:2);

end
end
end

function test = isPointInEllipse( eData, pData )

test = (((pData(1,1) - eData(1,1))^2 / eData(1,3)^2) + ((pData(1,2) -

eData(1,2))^2 / eData(1,4)^2)) <= 1;

end

function loadDataFileInMatrix( path, input, delimiterIn, headerlinesIn )

if nargin < 4
headerlinesIn = 1;
end

if nargin < 3
delimiterIn = ' ';
end

if ischar(input)

datastruct = importdata([path input], delimiterIn, headerlinesIn);

assignin('base',['dataset' regexprep(input,'\.(.)*','')],
datastruct.data(:,:));

elseif iscell(input)

for f=1:1:size(input,2)

datastruct = importdata([path input{f}], delimiterIn,

headerlinesIn);

assignin('base',['dataset' regexprep(input{f},'\.(.)*','')],
datastruct.data(:,:));

end

end

end

82
function [path, total_length] = opt2(path, distances)

n = size(path,2);

for i=1:n
for j=1:n-3
if change_in_path_length(path, i, j, distances) < 0
path = change_path(path, i, j);
end
end
end

total_length = distancePath(distances, path);

end

function [ optimalPath, optimalPathLength, optimalPathOrder, resultDist,

resultTime ] = optimization( fullloops, data, moveOptionsDPSO, moveOptionsPSO
)

resultDist = zeros(fullloops, 2);

resultTime = zeros(fullloops, 3);

optimalPathOrder = zeros(fullloops, size(data,1));

optimalPath = zeros(size(data,1), 2, fullloops);
optimalPathLength = Inf(fullloops,1);

parfor l=1:1:fullloops

tDpso = tic;
fprintf('DPSO\n');
[ path, total_length_dpso ] = psoOptDisc( data , 50,
moveOptionsDPSO);

dpsoResultTime = toc(tDpso) * 1000;

tPso = tic;
fprintf('PSO\n');
[ ~, total_length_pso, travelPoints ] = psoOpt( data, path, 50,
moveOptionsPSO);

psoResultTime = toc(tPso) * 1000;

resultDist(l,:) = [total_length_dpso total_length_pso];

resultTime(l,:) = [dpsoResultTime psoResultTime (dpsoResultTime +

psoResultTime)];

optimalPathOrder(l,:) = path;
optimalPath(:,:,l) = travelPoints;
optimalPathLength(l,1) = total_length_pso;
end
end

83
function plot_2d_function(func, lb, ub)
[X,Y] = meshgrid(lb(1):.1:ub(1), lb(2):.1:ub(2));
for (i=1:1:length(X(:,1)))
for (j=1:1:length(X(:,1)))
Z(i,j) = func([X(i,j); Y(i,j)]);
end
end
surfc(X,Y,Z)
end

function plotEllipse( x0, y0, r1, r2, color)

t = -pi:0.01:pi;

x = x0 + r1 * cos(t);
y = y0 + r2 * sin(t);

scatter(x0, y0, [], 'black', 'x', 'MarkerFaceColor', 'black')

if nargin < 5
plot(x, y, 'LineWidth', 2)
else
plot(x, y, 'LineWidth', 1, 'Color', color)
end

end

l = 95;

data = dataset_tspn2DE12_1;

travelPoints = optimalPath(:,:,l);

path = optimalPathOrder(l,:);

plotSpaceAfterPSO( data, path, travelPoints, true);

function plotSpace( data, path, plotPath )

if nargin < 3
plotPath = false;
end

figure;
hold all;

for i=1:1:size(data,1)
plotEllipse( data(i,1), data(i,2), data(i,3), data(i,4), 'k');
end

if plotPath == true

84
for j=1:1:size(path,2)-1
plot(data(path(1,j),1), data(path(1,j),2), 'o', 'MarkerSize', 5,
'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot(data(path(1,j+1),1), data(path(1,j+1),2), 'o', 'MarkerSize',
5, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot([data(path(1,j),1) data(path(1,j+1),1)],[data(path(1,j),2)
data(path(1,j+1),2)], 'b', 'LineWidth', 1.0)
end

plot([data(path(1,size(path,2)),1)
data(path(1,1),1)],[data(path(1,size(path,2)),2) data(path(1,1),2)], 'b',
'LineWidth', 1.0)
end

hold off
end

function plotSpaceAfterPSO( data, path, travelPoints, plotPath )

if nargin < 4
plotPath = false;
end

figure;
hold all;
for i=1:1:size(data,1)
plotEllipse( data(i,1), data(i,2), data(i,3), data(i,4) , 'k');
end

if plotPath == true
for j=1:1:size(path,2)-1
plot(travelPoints(path(1,j),1), travelPoints(path(1,j),2), 'o',
'MarkerSize', 5, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot(travelPoints(path(1,j+1),1), travelPoints(path(1,j+1),2),
'o', 'MarkerSize', 5, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot([travelPoints(path(1,j),1)
travelPoints(path(1,j+1),1)],[travelPoints(path(1,j),2)
travelPoints(path(1,j+1),2)], 'b', 'LineWidth', 1.0)
hold on
end

plot(travelPoints(path(1,size(path,2)),1),
travelPoints(path(1,size(path,2)),2), 'o', 'MarkerSize', 5,
'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot(travelPoints(path(1,1),1), travelPoints(path(1,1),2), 'o',
'MarkerSize', 5, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b');
plot([travelPoints(path(1,size(path,2)),1)
travelPoints(path(1,1),1)],[travelPoints(path(1,size(path,2)),2)
travelPoints(path(1,1),2)], 'b', 'LineWidth', 1.0)
end

hold off

end

85
function pRand = psoDiscRandomFact( particlePos, distances, iterSteps, type )

switch(type)
case'edgeExch'
pRand = edgeExchRandomFact( particlePos, iterSteps );

case'2opt'
pRand = twoOptRandomFact( particlePos, distances );

otherwise
pRand = particlePos;
end

end

function pRand = edgeExchRandomFact( particlePos, iterSteps )

if iterSteps <= 0 || isempty(iterSteps)

iterSteps = 2;
end

pRand = particlePos;
edgeExch = [];

for randIter=1:1:iterSteps
while true
startR = max(round((size(pRand,2)-1)*rand(1) + 1),1);
endR = startR + max(round((size(pRand,2)-1)*rand(1) + 1),1);

if r(startR, pRand) ~= r(endR, pRand)

edgeExch = [edgeExch; r(startR, pRand) r(endR, pRand)];
break
end
end
end
pRand = applyEdgeExchange( edgeExch, pRand );

end

function pRand = twoOptRandomFact( particlePos, distances )

[pRand, ~] = opt2(particlePos, distances);

end

func = @(x) ackley(x);

plot_2d_function(func, -5 * ones(2,1), 5 * ones(2,1));

function [ path, total_length, travelPoints ] = psoOpt( data, path,

swarmQuantity, moveOptions)

86
 anzCity = size(data,1);

particlePos = initSwarmMemberPso( data, swarmQuantity, 'random' );

personalBest = particlePos;

[ globalBest, ~ ] = findGlobalBestPso( path, particlePos );

lastVelocity = zeros(anzCity, 2, swarmQuantity);

if isfield(moveOptions,'initVelocity')
lastVelocity(:) = moveOptions.initVelocity;
else
lastVelocity(:) = 0.5;
end

noChangeCount = 0;

if isfield(moveOptions,'omega')
w = moveOptions.omega;
else
w = 0.5;
end

if isfield(moveOptions,'c1')
c1 = moveOptions.c1;
else
c1 = 0.8;
end

if isfield(moveOptions,'c2')
c2 = moveOptions.c2;
else
c2 = 0.5;
end

pi = 1;
while true

for p=1:1:size(particlePos,3)
for n=1:1:length(path)

phi1 = (1-0.1)*rand(1) + 0.1;

phi2 = (1-0.1)*rand(1) + 0.1;

newVelocity = w * lastVelocity(n,:,p) + phi1 * c1 *

(globalBest(n,:) - particlePos(n,:,p)) + phi2 * c2 * (personalBest(n,:,p) -
particlePos(n,:,p));

if isPointInEllipse( data(n,:), (particlePos(n,:,p) + newVelocity) ) == true

particlePos(n,:,p) = particlePos(n,:,p) + newVelocity;
lastVelocity(n,:,p) = newVelocity;
else

87
 XYproj = getNextBoarderPoint( data(n,:),
particlePos(n,:,p), (particlePos(n,:,p) + newVelocity) );

if ~isempty(XYproj)
deltaX = XYproj(1,1) - particlePos(n,1,p);
deltaY = XYproj(1,2) - particlePos(n,2,p);

if (deltaX / deltaX) > 0

signX = -1;
else
signX = 1;
end

if (deltaY / deltaY) > 0

signY = -1;
else
signY = 1;
end

lastVelocity(n,1,p) = 0;
lastVelocity(n,2,p) = 0;

particlePos(n,1,p) = XYproj(1,1) + (deltaX * signX) *

moveOptions.boundaryhandlingPercentage;
particlePos(n,2,p) = XYproj(1,2) + (deltaY * signY) *
moveOptions.boundaryhandlingPercentage;
end

end

distpersonalBest = distancePath( distance(personalBest(:,:,p)), path );

distAktPos = distancePath( distance(particlePos(:,:,p)), path
);
if distAktPos <= distpersonalBest
personalBest(n,:,p) = particlePos(n,:,p);
end
end
end

lastGlobalBest = globalBest;
[ globalBest, ~ ] = findGlobalBestPso( path, particlePos, globalBest );

if lastGlobalBest == globalBest
noChangeCount = noChangeCount + 1;
end

if isfield(moveOptions,'noChangeCountTh') && moveOptions.noChangeCountTh > 0

&& noChangeCount > moveOptions.noChangeCountTh
tfp = round( (swarmQuantity - 1) * rand(1) + 1 );
particlePos(:,:,tfp) = initSwarmMemberPso( data, 1 ,'random');
noChangeCount = 0;
end

88
if (isfield(moveOptions,'particleIter') && pi >= moveOptions.particleIter) ||
(isfield(moveOptions,'noChangeIterStop') && noChangeCount >
moveOptions.noChangeIterStop)
break
end

pi = pi + 1;
end
travelPoints = globalBest;

total_length = distancePath(distance(travelPoints), path);

end

function [ path, total_length ] = psoOptDisc( data, swarmQuantity,

moveOptions)

distances = distance(data(:,1:2));
particlePos = initSwarmMemberDpso( data, swarmQuantity, 'randomWithDup' );

personalBest = particlePos;

[ globalBestPath, ~ ] = findGlobalBestDpso( distances, particlePos );

noChangeCount = 0;

pi = 1;

while true

for p=1:1:size(particlePos, 1)

personalBest(p,:) = psoOptDiscMovement( particlePos(p,:),

personalBest(p,:), globalBestPath, distances, moveOptions );

distpersonalBest = distancePath( distances, personalBest(p,:) );

distAktPos = distancePath( distances, particlePos(p,:) );

if distAktPos < distpersonalBest

personalBest(p,:) = particlePos(p,:);
end
end

lastGlobalBest = globalBestPath;
[ globalBestPath, globalBestDist ] = findGlobalBestDpso( distances,
particlePos, globalBestPath );

if lastGlobalBest == globalBestPath
noChangeCount = noChangeCount + 1;
end

fprintf('%d\n', globalBestDist);

89
if (isfield(moveOptions,'particleIter') && pi >= moveOptions.particleIter) ||
(isfield(moveOptions, 'noChangeIterStop') && noChangeCount >
moveOptions.noChangeIterStop)
break
end

pi = pi + 1;
end
path = globalBestPath;
total_length = distancePath( distances, path );
end

function particlePos = psoOptDiscMovement( particlePos, personalBest,

globalBestPath, distances, moveOptions )

rLoc = (1-0)*rand(1) + 0;
rGlob = (1-0)*rand(1) + 0;
rRand = (1-0)*rand(1) + 0;
bRand = (1-0)*rand(1) + 0;

if moveOptions.bLoc >= 0 && moveOptions.bLoc <= 1

bLoc = moveOptions.bLoc;
else
bLoc = 0.3;
end

if moveOptions.bGlob >= 0 && moveOptions.bGlob <= 1

bGlob = moveOptions.bGlob;
else
bGlob = 0.3;
end

edgeExch = getListOfEdgeExchanges( particlePos, personalBest );

if ~isempty(edgeExch) && round(size(edgeExch,1) * rLoc * bLoc) > 0
dLoc_list = edgeExch(1:round(size(edgeExch,1) * rLoc * bLoc), :);
else
dLoc_list = [];
end
dLoc = applyEdgeExchange( dLoc_list, particlePos );
edgeExch = getListOfEdgeExchanges( particlePos, globalBestPath );
if ~isempty(edgeExch) && round(size(edgeExch,1) * rGlob * bGlob) > 0
dGlob_list = edgeExch(1:round(size(edgeExch,1) * rGlob * bGlob), :);
else
dGlob_list = [];
end
dGlob = applyEdgeExchange( dGlob_list, particlePos );
edgeExch = getListOfEdgeExchanges( dGlob, dLoc );
if ~isempty(edgeExch) && round(size(edgeExch,1) * 0.5) > 0
d_temp_list = edgeExch(1:round(size(edgeExch,1) * 0.5), :);
else
d_temp_list = [];
end
d_temp = applyEdgeExchange( d_temp_list, dGlob );

90
if strcmp(moveOptions.randArt, 'randomStart')
pRand = psoDiscRandomFact( particlePos, distances,
moveOptions.vRandIter, moveOptions.vRandType );
elseif strcmp(moveOptions.randArt, 'randomTemp')
pRand = psoDiscRandomFact( d_temp, distances, moveOptions.vRandIter,
moveOptions.vRandType );
else
pRand = psoDiscRandomFact( particlePos, distances,
moveOptions.vRandIter, moveOptions.vRandType );
end

edgeExch = getListOfEdgeExchanges( d_temp, pRand );

if ~isempty(edgeExch) && round(size(edgeExch,1) * rRand * bRand) > 0
v_rand_list = edgeExch(1:round(size(edgeExch,1) * rRand * bRand), :);
else
v_rand_list = [];
end

particlePos = applyEdgeExchange( v_rand_list, d_temp );

end

function result = r(idx, path)

pathLength = size(path, 2);

if idx > pathLength

result = idx - pathLength;
else
result = idx;
end
end

function total_length = total_length_of_cycle(distances, path)

total_length = 0;
for i = 1:1:length(path)-1
total_length = total_length + distances(path(i), path(i+1));
end
total_length = total_length + distances(path(end), path(1));
end

91