ruler scale

clear
close all
rgb = imread('Scan_Pic0003.jpg');
imshow(rgb);
Warning: Image is too big to fit on screen; displaying at 25% 
load pic_data
I2C = imcrop(rgb,RECT);
imshow(I2C);
Warning: Image is too big to fit on screen; displaying at 50% 

filter it

h = fspecial('gaussian',5,2);
I2 = imfilter(I2C,h);
imshow(I2)
Warning: Image is too big to fit on screen; displaying at 50% 
y = g(14,:)';
ye = zeros(4096,1);
offset = floor((4096-length(y))/2)
ye(offset+(1:length(y)))=y;
freq = -2048:2047;
yf = fft(ye);
yf = fftshift(abs(yf));
plot(freq,yf);
offset =

   572

idx = find(abs(freq-170)<10);
f = freq(idx)';
yfc = yf(idx);
[ys is] = sort(yfc,'descend');
idx = [-1:1] + is(1);
c = polyfit(f(idx),yfc(idx),2); % quadratic fit
hmax = -c(2)/(2*c(1)) % harmonics
dpmm = 4096/hmax
dpi =  dpmm*25.4
hmax =

  173.6420


dpmm =

   23.5888


dpi =

  599.1546

idx = [-2:1] + is(1);
fx = linspace(min(f(idx)),max(f(idx)),51);
yfit = polyval(c,fx);
plot(f,yfc,fx,yfit);
sz = size(I2);
y = im2double(I2(floor(sz(1)/2),:,1));
xc = floor(sz(2)/2);
idx = (-80:80);
plot(idx,y(xc+idx));
title(sprintf('center pixel %d\n',xc));
axis([-80 80 -0.1 1.1]);
xlabel('pixels');
ylabel('relative brightness');
dpmm = 600/25.4;
[ys is] = sort(y<0.4,'descend');
npts = floor(sz(2)/dpmm);
xo = floor(is(1));
xcen = zeros(1,npts);
ptn = zeros(1,npts);
x = linspace(-3,3,101);
N=0;
jdx = (-10:10);
idx = (-3:3);
for k=1:npts
    [ys is] = min(y(xo+jdx));
    xo = floor(xo+is(1)-11);
    c = polyfit(idx,y(xo+idx),2);
    dx = -c(2)/(2*c(1));
    yfit = polyval(c,x);
%     plot(jdx,y(xo+jdx),x,yfit);
%     title(sprintf('pt %d center %d',k,xo));
    if (abs(dx)<3)
        N = N + 1;
        ptn(N) = k;
        xcen(N) = xo + dx;
        xo = floor(xcen(N) + dpmm);
    else
        xo = floor(xo+dpmm);
    end
    if (xo+18)> sz(2)
        break;
    end
end
xcen = xcen(1:N);
ptn = ptn(1:N);
npts = N;
plot(y(1:200),'k');
idx = find(xcen<200);
hold on;
for k=idx
    plot([xcen(k) xcen(k)],[0 1]);
end
hold off
c = polyfit(ptn,xcen,1);
fprintf('c(1) %g c(2) %g\n',c(1:2));
xfit = polyval(c,ptn);
plot(ptn,xcen-xfit);
xlabel('point');
ylabel('error');
grid;

dpi = c(1)*25.4
c(1) 23.5987 c(2) -5.67507

dpi =

  599.4060