I am working on images to detect and fit the largest possible circle in any of the free areas of an image containing distributed particles: 

(able to detect the location of particle).

One direction is to define a circle touching any 3-point combination, checking if the circle is empty, then finding the largest circle among all empty circles. However, it leads to a huge number of combination i.e. C(n,3), where n is the total number of particles in the image.

I would appreciate if anyone can provide me any hint or alternate method that I can explore. 

Lets do some maths my friend, as maths will always get to the end!

Wikipedia:

In mathematics, a Voronoi diagram is a partitioning of a plane into regions based on distance to points in a specific subset of the plane. 

For example: 

rng(1)
x=rand(1,100)*5;
y=rand(1,100)*5;

voronoi(x,y);

The nice thing about this diagram is that if you notice, all the edges/vertices of those blue areas are all to equal distance to the points around them. Thus, if we know the location of the vertices, and compute the distances to the closest points, then we can choose the vertex with highest distance as our center of the circle.

Interestingly, the edges of a Voronoi regions are also defined as the circumcenters of the triangles generated by a Delaunay triangulation.

So if we compute the Delaunay triangulation of the area, and their circumcenters 

dt=delaunayTriangulation([x;y].');
cc=circumcenter(dt); %voronoi edges

And compute the distances between the circumcenters and any of the points that define each triangle: 

for ii=1:size(cc,1)
    if cc(ii,1)>0 && cc(ii,1)<5 && cc(ii,2)>0 && cc(ii,2)<5
    point=dt.Points(dt.ConnectivityList(ii,1),:); %the first one, or any other (they are the same distance)
    distance(ii)=sqrt((cc(ii,1)-point(1)).^2+(cc(ii,2)-point(2)).^2);
    end
end

Then we have the center (cc) and radius (distance) of all possible circles that have no point inside them. We just need the biggest one! 

<pre.[r,ind]=max(distance); %Tada!

Now lets plot 

hold on

ang=0:0.01:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);

point=cc(ind,:);

voronoi(x,y)
triplot(dt,'color','r','linestyle',':')
plot(point(1)+xp,point(2)+yp,'k');
plot(point(1),point(2),'g.','markersize',20);

Notice how the center of the circle is on one vertex of the Voronoi diagram. 

NOTE: this will find the center inside [0-5],[0-5]. you can easily modify it to change this constrain. You can also try to find the circle that fits on its entirety inside the interested area (as opposed to just the center). This would require a small addition in the end where the maximum is obtained.