Registration using particle filter in matlab

I wrote the following code as an implementation to the registration algorithm presented in an article named:"E. Arce-Santana, D. Campos-Delgado and A. Alba, Affine image registration guided by particle filter, IET Image Process. 6 (2012)".

Question 1: When I run the code no matter how many particles or iterations I choose, the output is still inaccurate (by inaccurate I mean even with increasing the number of particles that according to the algorithm should increase the performance, I can not reach that). I do not know what is the problem?

Question 2: I commented a formula for updating the particle's weight, my question is, does I implement the equation write or the one that appears above it (which relied on the entropy) is the right one and therefore delete the commented one and leave the entropy based equation?

The code is as follow:

%%clear everything
clear all
close all
clc

%%Read the image data

% I1: Reference Image
I1=imread('cameraman.tif');
I1=double(imresize(I1,[256 256]));
figure,imshow(I1)
[I1r I1c I1d]=size(I1);

%I2: Target image
I2=randomtransform(I1); %user-defined function to generate random transformation 
% I2=double(imresize(I2,[I1r I1c]));
figure,imshow(I2)

%%Particle Filter Steps:
%%%%%Input:
n=4;         % Related to the initialization
m=256;       % Related to the initialization
N=20;       %(no. of iteration)
M=100;       %(no. of particles)
phai=[];     %(state vector of the affine transformation parameters)
w=[];        %(the vector of the phai weights)
%k: iteration no.
%i: particle no.
Vk=1;        % noise variance in the system (Used in the predection state)
w_v=1;       % weights update equation related variance
beta=0.99;   % annealing factor
%%%%%Output
phai_est=[]; %(Final estimated state vector of affine parameters)

%%%%%Steps:
%%%Step 1: Random generation of the particles
% The ranges are obtained from the paper
    rotationAngle=double(int8(unifrnd(-pi/4,pi/4,1,1)));
    scalingCoefX = double(unifrnd(0.5,2,1,1));
    scalingCoefY = double(unifrnd(0.5,2,1,1));
    shearingCoefX = double(unifrnd(0.5,2,1,1));
    shearingCoefY = double(unifrnd(0.5,2,1,1));
    translationCoefX = double(int8(unifrnd(-I1r/2,I1r/2,1,1)));
    translationCoefY = double(int8(unifrnd(-I1c/2,I1c/2,1,1)));
    %The initialization of the first phai
    phai(1,:)=[round(rotationAngle*10^n)/10^n, round(scalingCoefX*10^n)/10^n, round(scalingCoefY*10^n)/10^n, round(shearingCoefX*10^n)/10^n, round(shearingCoefY*10^n)/10^n, round(translationCoefX*10^n)/10^n, round(translationCoefY*10^n)/10^n]';
    %Make the randomly generated particles from the initial prior gaussian distribution
    for i = 1:M
    phai(i,:) = phai(1,:) + sqrt(2) * randn; %2: the variance of the initial esimate
    w(i,:)=1/M;
    end
    % Normalize the weights:
    w = w./sum(w);

%%%Step2: Resample process
for k=1:N
   for i=1:M
       % rand: u (Uniform random value between 0 and 1
       j=find((cumsum(w) >= max(w)),1,'first');
       phai_select(i,:)=phai(j,:);
       phai(i,:)=phai_select(i,:)+(sqrt(Vk^2)*randn);
       I2_new=targettransform(I2,phai(i,:)); %user-defined function to apply the generated transformation to the target image 
       E_I1=entropy(I1);
       I=E_I1+entropy(I2_new)-joint_entropy(I1,I2_new); %joint_entropy: user defined function to calculate joint entropy of the two images
       w(i)=(1/sqrt(2*pi*w_v))*exp(-((E_I1-I)^2)/(2*w_v));
%        w(i)=prod(prod(((1/sqrt(2*pi*w_v))*exp(-((I2_new-I1)^2)/(2*w_v)))));
   end
       % Normalize the weights
       w = w./sum(w);
       % Reduce the noise standard deviation
       Vk=beta*Vk;
       phai_est=mean(phai);
end