I am a masters student working on a project for one of my professors and i am coming across an issue. I am working on a 6 variable GA scheme for Digital image correlation. I have created a majority of the code but i want to be able to add in the ability to to subsize imaging. the code right now will only return hole numbers for the amount of displacement in the u and v direction. i would like to be able to have the ability to have fractions of displacement in terms of the pixels. i have the code that will interpolate the image to make the number of pixels smaller but i am having trouble inserting it into the code, keep coming up with errors. i would like someone who is more knowledgable than me with matlab to see if they could help. i know that the interpolation needs to go inside the while loop when it is looking at the the subset size. if anyone could help please let me know. here is the code for the DIC
clear; clc;
u_displacement = 5;
v_displacement = 6;
size = 100;
%%%%%%% generate artificial array a1 and a2 %%%%%%%%%%
IMAGE = imread('8a.bmp');
a1 = double(IMAGE(1:size,1:size));
for i = 1:size
for j = 1:size
newrow(i,j) = i + v_displacement + round((j-1)*0.1 + (i-1)*0.2);
newcol(i,j) = j + u_displacement + round((j-1)*0.3 + (i-1)*0.4);
end
end
for i = 1:size
for j = 1:size
a2(newrow(i,j),newcol(i,j)) = a1(i,j);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% a1 and a2 generated %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tpointrow = 1;
tpointcol = 1;
num_lookup_table = 2000; % number of values for u and v generated in the look up tables
subset_size = 60;
ipop = 500; % initial population of u and v indices
iguesu_lower = 0; % initial guess for lower limit of u
iguesu_upper = 15; % initial guess for upper limit of u
iguesv_lower = 0; % initial guess for lower limit of v
iguesv_upper = 10; % initial guess for upper limit of u
generations = 100;
crossovertimes = 3;
crossoverrate = 0.5;
wpop = 100;
b = 0.5;
mutationrate = 0.5;
fact = subset_size;%ceil(subset_size/2); %%%factor to be added and subtracted to point of interest
%%%% generate initial table of values for u and v %%%%%%%
for i = 1:(num_lookup_table+1)
ulookup1(i) = floor(rand*(iguesu_upper-1.5))+iguesu_lower;
vlookup1(i) = floor(rand*(iguesv_upper-1.5))+iguesv_lower;
s11(i) = (round(rand*0.5*10))/10;
s21(i) = (round(rand*0.5*10))/10;
s31(i) = (round(rand*0.5*10))/10;
s41(i) = (round(rand*0.5*10))/10;
end
[ulookup indx_ulookup] = sort(ulookup1);
[vlookup indx_vlookup] = sort(vlookup1);
[s1lookup indx_s1lookup] = sort(s11);
[s2lookup indx_s2lookup] = sort(s21);
[s3lookup indx_s3lookup] = sort(s31);
[s4lookup indx_s4lookup] = sort(s41);
for k = 1:ipop
u(k) = ceil(rand*num_lookup_table)+1;
v(k) = ceil(rand*num_lookup_table)+1;
ps1(k) = ceil(rand*num_lookup_table)+1;
ps2(k) = ceil(rand*num_lookup_table)+1;
ps3(k) = ceil(rand*num_lookup_table)+1;
ps4(k) = ceil(rand*num_lookup_table)+1;
end
%%%generate initial pop%%%%%%%%%%
for k = 1:ipop
num = 0;
denum1 = 0;
denum2 = 0;
corr1(k)= 0;
for row = (tpointrow):(tpointrow+fact)
for col = (tpointcol):(tpointcol+fact)
DX = (col - tpointcol);
DY = (row - tpointrow);
dudx = s1lookup(ps1(k));
dudy = s2lookup(ps2(k));
dvdx = s3lookup(ps3(k));
dvdy = s4lookup(ps4(k));clc;
r = row+vlookup(v(k))+round(dudx*DX+dudy*DY);
c = col+ulookup(u(k))+round(dvdx*DX+dvdy*DY);
num = num + a1(row,col)*a2(r,c);
% a1(row,col)
% a2(r,c)
denum1 = denum1 + (a1(row,col))^2;
denum2 = denum2 + (a2(r,c))^2;
end
end
corr1(k) = -num/sqrt(double(denum1)*double(denum2));
end
[sortedcost1 indx] = sort(corr1);
newu = u;
newv = v;
news1 = ps1;
news2 = ps2;
news3 = ps3;
news4 = ps4;
g = 1;
flag = 0;
corr = sortedcost1(1:wpop);
bestu(1) = newu(indx(1));
bestv(1) = newv(indx(1));
bests1(1) = news1(indx(1));
bests2(1) = news2(indx(1));
bests3(1) = news3(indx(1));
bests4(1) = news4(indx(1));
trackcost(1) = sortedcost1(1);
if (corr(1) <=-0.99)
flag = 1;
end
while (g<=generations)&(flag~=1)
newu1 = newu; clear newu;
newv1 = newv; clear newv;
news11 = news1; clear news1;
news21 = news2; clear news2;
news31 = news3; clear news3;
news41 = news4; clear news4;
for i = 1:(wpop)
newu(i) = newu1(indx(i));
newv(i) = newv1(indx(i));
news1(i) = news11(indx(i));
news2(i) = news21(indx(i));
news3(i) = news31(indx(i));
news4(i) = news41(indx(i));
end
clear sortedcost indx;
for crosstimes = 1:crossovertimes
if(rand<=crossoverrate)
indextocrossoveru = ceil(rand*wpop);
indextocrossoverv = ceil(rand*wpop);
indextocrossovers1 = ceil(rand*wpop);
indextocrossovers2 = ceil(rand*wpop);
indextocrossovers3 = ceil(rand*wpop);
indextocrossovers4 = ceil(rand*wpop);
swap = newu(indextocrossoveru);
newu(indextocrossoveru) = newv(indextocrossoverv);
newv(indextocrossoverv) = swap;
swap = news1(indextocrossovers1);
news1(indextocrossovers1) = news2(indextocrossovers2);
news2(indextocrossovers2) = swap;
swap = news3(indextocrossovers3);
news3(indextocrossovers3) = news4(indextocrossovers4);
news4(indextocrossovers4) = swap;
swap = news2(indextocrossovers2);
news2(indextocrossovers2) = news4(indextocrossovers4);
news4(indextocrossovers4) = swap;
end
end
for i = ((wpop/2)+1):2:(wpop-1)
newu(i) = ceil((b*newu(i-(wpop/2))+(1-b)*newu(i-(wpop/2-1))));
newu(i+1) = ceil(((1-b)*newu(i-(wpop/2))+(b)*newu(i-(wpop/2-1))));
newv(i) = ceil((b*newv(i-(wpop/2))+(1-b)*newv(i-(wpop/2-1))));
newv(i+1) = ceil(((1-b)*newv(i-(wpop/2))+(b)*newv(i-(wpop/2-1))));
news1(i) = ceil((b*news1(i-(wpop/2))+(1-b)*news1(i-(wpop/2-1))));
news1(i+1) = ceil(((1-b)*news1(i-(wpop/2))+(b)*news1(i-(wpop/2-1))));
news2(i) = ceil((b*news2(i-(wpop/2))+(1-b)*news2(i-(wpop/2-1))));
news2(i+1) = ceil(((1-b)*news2(i-(wpop/2))+(b)*news2(i-(wpop/2-1))));
news3(i) = ceil((b*news3(i-(wpop/2))+(1-b)*news3(i-(wpop/2-1))));
news3(i+1) = ceil(((1-b)*news3(i-(wpop/2))+(b)*news3(i-(wpop/2-1))));
news4(i) = ceil((b*news4(i-(wpop/2))+(1-b)*news4(i-(wpop/2-1))));
news4(i+1) = ceil(((1-b)*news4(i-(wpop/2))+(b)*news4(i-(wpop/2-1))));
end
%%%%% mutate%%%%%%
if(rand<=mutationrate)
newu(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
newv(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
news1(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
news2(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
news3(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
news4(ceil(rand*wpop)) = floor(rand*(num_lookup_table-1.5))+1;
end
%%%%%%%%%%
%%%%%%%% objective function
num = 0;
denum1 = 0;
denum2 = 0;
% if(rand>0.5)
% sign = -1;
% else
sign = 1;
% end
for k = 1:wpop
corr(k) = 0;
for row = (tpointrow):(tpointrow+fact)
for col = (tpointcol):(tpointcol+fact)
DX = (col - tpointcol);
DY = (row - tpointrow);
dudx = s1lookup(news1(k));
dudy = s2lookup(news2(k));
dvdx = s3lookup(news3(k));
dvdy = s4lookup(news4(k));
r = round(row+vlookup(newv(k))+dudx*DX+dudy*DY);
c = round(col+ulookup(newu(k))+dvdx*DX+dvdy*DY);
num = num + a1(row,col)*a2(r,c);
denum1 = denum1 + (a1(row,col))^2;
denum2 = denum2 + (a2(r,c))^2;
end
end
corr(k) = -num/sqrt((denum1)*(denum2));
end
[sortedcost indx] = sort(corr);
trackcost(g+1) = sortedcost(1);
bestu(g+1) = newu(indx(1));
bestv(g+1) = newv(indx(1));
bests1(g+1) = news1(indx(1));
bests2(g+1) = news2(indx(1));
bests3(g+1) = news3(indx(1));
bests4(g+1) = news4(indx(1));
if (sortedcost(1) <=-0.99)
flag = 1;
end
g = g+1;
end
if (flag == 1)
fprintf('Correct displacements were predicted after %d generations\n',(g))
fprintf(' u = %d pixels v = %d pixels\ns1 = %5.3f s2 = %5.3f s3 = %5.3f s4 = %5.3f',ulookup(bestu(g)),vlookup(bestv(g)),s1lookup(bests1(g)),s2lookup(bests2(g)),s3lookup(bests3(g)),s4lookup(bests4(g)))
else
fprintf('Correct displacements could not be predicted with sufficient accuracy even after %d generations\n',(g-1))
end
figure (1); plot(trackcost,' - o');
title('movement of correlation');
for count = 1:numel(bestu)
toplotu(count) = ulookup(bestu(count));
toplotv(count) = vlookup(bestv(count));
toplots1(count) = s1lookup(bests1(count));
toplots2(count) = s2lookup(bests2(count));
toplots3(count) = s3lookup(bests3(count));
toplots4(count) = s4lookup(bests4(count));
end
figure (4); plot(toplotu,' - *');
title('movement of x-displacement');
figure (5); plot(toplotv,' - *');
title('movement of y-displacement');
figure (6); plot(toplots1,' - o');
title('movement of s1');
figure (7); plot(toplots2,' - o');
title('movement of s2');
figure (8); plot(toplots3,' - o');
title('movement of s3');
figure (9); plot(toplots4,' - o');
title('movement of s4');
figure(10)
imshow(uint8(a1));
figure(11);
imshow(uint8(a2));
and here is the code for the interpolation (it does alot more than i need but the part of it that does the interpolation to create, lets say a 20 x 20 into a 200 x 200 pixel, is what i need.
% COR.M
% Digital Image Correlation File
% % Uses a tif file - 'speckle1.tif'
%
clear;
%
speckle=imread('speckle1.tif');
mesh(speckle);
pause;
%imshow(speckle);pause;
%imhist(speckle);pause;
s=size(speckle)
% change int to double to work with mesh
% to change back use uint8
a=double(speckle);
s1=size(a)
mesh(a),rotate3d on;
pause;
%surf(a), %b = a > 128;
%imshow(b); pause;
c=1000+a;
mesh(c),rotate3d on;
axis([0 300 0 300 0 1400]);
pause;
% Looking at a 10x10 area
d=zeros(240,320);
roi=zeros(10,10);
start_x=100;
start_y=100;
for i=start_x:start_x+9
for j=start_y:start_y+9
d(i,j)=1500;
roi(i-start_x+1,j-start_y+1)=c(i,j);
end
end
mesh(roi);
axis([1 10 1 10 0 1400]);
pause;
e=c+d;
mesh(e),rotate3d on;
axis([0 300 0 300 0 3000]);
pause;
%
% This increases resolution of the region of interest (100x100)
f=zeros(100,100);
for i=0:10:90
for l=0:10:90
for j=1:10
for k=1:10
f(j+i,k+l)=c(start_x+i/10,start_y+l/10);
end
end
end
end
mesh(f);pause;
%
%Using Interpolation (Can use 'bilinear', 'bicubic' or 'nearest')
%
xa=1:1:10;
ya=1:1:10;
[x,y]=meshgrid(xa,ya);
z=roi;
mesh(z);pause;
xa1=1:0.1:10;
ya1=1:0.1:10;
[x1,y1]=meshgrid(xa1,ya1);
z1=interp2(x,y,z,x1,y1,'bilinear'); %default is 'bilinear'
mesh(z1);
pause;