%%

% 
% # 1.Data Sources dataset has been downloaded. Tiny Images dataset, which 
%consists of 79,302,017 images, each being a 32x32 color image. This data is
%stored in the form of large binary files which can be accesed by a Matlab toolbox 
%that we have written. (i) the images themselves; (ii) their associated metadata 
%(filename, search engine used, ranking etc.);   
load('tiny_index.mat');
% The user marks an ellipse at approximate size and location of the object
% with mouse click
im=imread('images.png')
figure, imshow(im);
h = imellipse(gca);
addNewPositionCallback(h,@(p) title(mat2str(p,3)));
fcn = makeConstrainToRectFcn('imellipse',get(gca,'XLim'),get(gca,'YLim'));
setPositionConstraintFcn(h,fcn);
color = getColor(h);
object = getPosition(h);
%objcolor=getColor(object)
obj= object([1 2]);
setColor(h,'yellow')
% x is identified by clicking mouse and then clicking return. x is part of object or part of background
% color of x
[x, y] = getpts(gca);
x=[x,y];
%computing color histogram for original image
new=rgb2gray(im);
figure
imshow(new);
imhist(new)
%%%%%%%%%%%%%%
I = imread('images.png');
imshow(I)
%threshold the image in order to produce a binary input image
bw = im2bw(I, graythresh(I));
imshow(bw)
%Fill in the holes.
bw2 = imfill(bw, 'holes');
imshow(bw2)
%You can easily see that there are 10 objects in the image. But to start
%doing geometric mefminconasurements, we have to be able to identify which foreground 
%pixels belong to which object. 
L = bwlabel(bw2);
%L is a label matrix. It's the same size as bw2. It contains the value 0 in each
%location that's a background pixel. It contains positive integer values in locations
%corresponding to a labeled foreground object. 
imshow(L == 23)
title('Object 23')
Inp=rgb2gray(I);
subplot(2,1,1);
imhist(Inp)
title('originalimage''s color histogram')
grayImg=mat2gray(L);
subplot(2,1,2);
imhist(grayImg)
title('object''scolor histogram')

%f(x) denotes the log likelihood ratio for observing the color I(x) at
%apoint x
f(x)=log(imhist(Inp))-log(imhist(grayImg));
i=1; beta=0.5;
 %s is a gaussian smoothed version of input image
s(x)= imgaussfilt(f(x),i);
%g is measure of local edge strength or photo consistency in multiview
%reconstruction with suitable choices of parameters a and sigma
g(x)=1./(1+beta.*abs(gradient(s(x))).^2);
%pixelposition=getpixelPosition(I)
%S is hyper surface (set of closed boundaries or closed surfaces.E(S) is sum of integrals of f and g that are defined informally 
%as the signed area of the region in the xy-plane that is bounded by the
%graph of f and g, the x-axis and the vertical lines x = int(S) and x = S respectively
F=integral(f(x),int(S),0);
%the integral of the Gaussian curvature g(x) over the entire surface S is determined by 
%the Euler characteristic
G=2*3.14159.*X(S); 
E(S)=sum(F,G);
%E(S) functionals can be globally optimized by means of convex relaxation and thresholding
%u(x) is indicator function representing hyper surface S on binary functions of bounded variation
% it is interior of S if it equals to 1 and it is exterior of S if it equals
% to 0.find u that minimizes the function E(u) over the set bounded
% variation BV[0,1]
if u(x)<0
    u(x)=0;
else if u(x)>=0 && u(x)<=1
        u(x)=u(x);
    else
        u(x)=1;
    end
end
cvx_begin 
% minimizing energy E(u)
options = optimoptions('fminunc','Algorithm','trust-region','GradObj','on');  
%Find the location and objective function value of the minimizer u*
%starting at u0=[-1,2]
[u, ~] = fmincon(@rosenbrock,[-1,2],options);
cvx_end
%Euler_lagrange equation can be solved. 
diff(g.*(gradient(u)./abs(gradient(u))))-f-0;
%gradient is less than the default value of the function tolerance.
lambda=[1;0.5;0.3;1.5;0.6;1;1.8;2;0.2;0.16;1.4;2];
%the moments on segementation could be imposed by area of the shape u to be
%bounded by a constant
cvx_begin
variable u(x)
%c is constant constraining the volume of shape in 3 dimensional domain
options = optimoptions('fmincon','Algorithm','interior-point','MaxIter',1500); 
[x,~] = fmincon(@areaConstraint,[0,0],[],[],[],[],[0,0],[1,0],@unitdisk,options);
cvx_end
%Euler-lagrange equation is computed by derivation of E0(u,lambda_0) for
%the area constraint. constraint parameter c is divided by the size of the
%image domain integral over omega equals to absolute value of omega in
%order to distribute the update value over all points
diff(g.*(gradient(u)./abs(gradient(u)))-f-2.*lambda(1).*(int(u)-c./integral()))=0;
% in the case of moment hard constraints, all constraints are enforced 
% during the optimization by back-projecting the current segmentation onto
% the constraint set. segmentation shape u are projected onto the
% intersection of two convex sets: convex set[0,1] and intersection of
% moment constraints
cvx_begin
variable u0(n)% for any u0, projection u onto area constraint has to solve
%the following convex program
variable v(n)
u=argmin (1./2).*norm(v-u0,2).^2;
subject to
         le(v(x)-1,0)
         le(-v(x),0)
         norm(v,1)-c-0;
cvx_end
%The KKT conditions are analogous to the condition that the gradient must be 
%zero at a minimum, modified to take constraints into account. The KKT conditions
%use the auxiliary Lagrangian function:.
%%
% 
% $$L(x,?)=f(x)+?.*g(x)+?.*h(x)$$
% 
diff(g.*(gradient(u)./abs(gradient(u)))-f-2.*lambda(1).*(int(u)-c./integral()))=0;
u=u0;
v=0;%opyimal value needed to be added to u to fulfill area constraint 
%%
% 
% $$int(u)=c$$
% 
%and u belongs to [0,1]
g0=0;
g1=0;
Dc=c-norm(u0);
x1=min(u0(x)<1);
x2=min(u0(x)<0);
while Dc>0 
    D1=1-u0(x1)+v;
    D2=-u0(x2)+v;
    if D1>D2
        v=v-min(D2,Dc./x2-x1);
        Dc=Dc-min((x2-x1).*D2,Dc);
        x2=min((u0(x)-v)<1);
    else
v=v-min(D2,Dc./x2-x1);
Dc=Dc-min((x2-x1).*D1,Dc);
x1=min((u0(x)-v)<1);
    end
end
 %The KKT conditions are:  
 u(x)-u0(x)+g1(x)-g0(x)+v-0;
 intersect(le(u(x),1),ge(g1(x),0))
intersect(ge(u(x),0),ge(g0(x),0))
union(g1(x)-0,u(x)-1);
union(g0(x)-0, u(x)-0);
norm(u,1)=c;
% imposing the centroid of shape as a soft constraint by minimizing energy E1
cvx_begin
variable u(x)
variable m(n)% m is centre 
m=integral(x.*u)./integral(u);
options = optimoptions('fmincon','Algorithm','interior-point','MaxIter',1500); 
[x,~] = fmincon(@centroidConstraint,[0.25,0.25],[],[],[],[],[0,0],[1,0],@constraintFun,options);
cvx_end
%the centroid soft constraint yields an Euler-lagrange equation 
for i=2:d
    % * m(i) and x(i)denote i-th elements of vectors m and x.
diff(g.*(gradient(u)./abs(gradient(u)))-f-2.*lambda(i).*int((m(i)-x(i)).*u.*(m(i)-x(i))))=0;
end
%segmentation shape u are projected to the nearest u* in set
%%
% 
% $$l1= int((m(i)-x(i)).*u $$
% 
cvx_begin
s=argmin(norm(z,1));
subject to 
for i=1:26
 l1=int((m(i)-x(i)).*u)-0;
M=int((m(i)-x(i)).*z)+ l1;
end
cvx_end
% * the corresponding soft constraint on the covariance matrix can be derived 
%by adding a respective term to the original energy and imposed by 
% * considering the following convex set 
%%
% 
% $$integral(norm(x-m,2).*u)./integral(u)=A$$
% 
cvx_begin
variable x(n)
variable u(n)
variable m(n)% m is centre that is intersection of A and 
%%
% 
% $$m=integral(x.*u)./integral(u)$$
% 
variable A(i,j)% a symmetric matrix
A=integral(norm(x-m,2).*u)./integral(u);
m=integral(x.*u)./integral(u);
options = optimoptions('fmincon','Algorithm','interior-point','MaxIter',1500); 
[x,fval] = fmincon(@covarianceobjFun,[0.25,0.25],[],[],[],[],[0,0],[1,0],@covarianceconFun,options);
cvx_end
%%
% # 2. problem formulation: The problem of segmentation of a given image is considered.
%Although it has been demonstrated that imposing shape priors improve image
%segmentation,  many existing approaches have some problems, for example,
%resulted solutions are locally optimal.Also, in medical image
%segmentation, malformation og organs must be detected .unfortunately,deviations of observd
%from training shapes are suppressed by  shape prior 
% 
% # 3.proposed solution: to solve this problem,a new approach has been
% suggested that user imposes constraints (area, centroid, and
% covariance) as convex ones on moments of arbitrary order within optimization
%framework but without imposing on local shape . image segmentations can be
%computed by minimizing convex functionals E(u) over convex set bonded variation
%[0,1]under respective convex constraints then, Euler-lagrange equations can
%be computed using gradient descent
%%
%  # 4.optimization tools that are used to solve this optimization problem are 
%Unconstrained and Constrained nonlinear optimization problems. Unconstrained
%nonlinear composed of only objective functions while Constrained nonlinear
%are composed of linear or nonlinear objective functions and may be subject
%to linear and nonlinear constraints. The interior point algorithm was used 
%to solve minimization problem.
 % 
%%
% # 5. Conclusion and analysis: In conclusion, occurance of Inf's and NaN's
% was frequent, As a result, Ihave been able to reproduce some of results reported in the original paper.
% through writing the algorithm, I have tried make the algorithm beave as
% described by authers. one of drawbacks of proposed solution is that I
% couldnot connect dataset and minimization problems
% 
% 
% # 6. list of sorce files and functions necessary for implementation:
% many functions have been created to implement this algorithm :
% rosenbrock,areaConstraint,unitdisk,centroidConstraint,constraintFun,
%CovarianceobjFun, and covarianceconFun. In addition, many available
%functions are used, such as optimoptions, fmincon,fminunc, diff, integral,
%imshow, imhist, rgb2gray,mat2gray and many others http://www.mathworks.com

% 
