addpath /research/src/DNAdesign
addpath /research/src/AFMutils

cd ~winfree/Research/DNA/AFM_old/


figure(1); [im,in]=showafm('data.rizal/02-07-20/QBOX.097',0); im=flipud(im);
figure(1); imagesc(im); axis square

figure(2)
mask = im(290:316,185:199);
mask = (mask+fliplr(mask)+flipud(mask)+fliplr(flipud(mask)))/4;
imagesc(mask); colormap gray; axis square

figure(3)
imc = conv2(im,mask,'same'); imagesc(imc); colormap gray; axis square

figure(2)
mask=imc(218:244,208:222);  
mask = (mask+fliplr(mask)+flipud(mask)+fliplr(flipud(mask)))/4;
imagesc(mask); colormap gray; axis square   
figure(3)
imc = conv2(im,mask,'same'); imagesc(imc); colormap gray; axis square
%imagesc(imc>1000); 

%% or we could smooth the (already smooth) mask...
figure(2)
x=-10:10; g=(exp(-x.^2/4)/2); plot(x,g) 
mask2=conv2(g,[1],mask,'same'); imagesc(mask2); axis square
figure(3)
imc = conv2(im,mask2,'same'); imagesc(imc); colormap gray; axis square
imc = imc/max(max(imc));
%imagesc(imc>.2);                                                 
%% does look a tad better.

%erode_imc= bwmorph(imc>.25, 'shrink', 10); figure(2); imagesc(erode_imc); colormap gray; axis square
%erode_imc= bwmorph(imc>.15, 'shrink', 5); figure(3); imagesc(erode_imc); colormap gray; axis square
%imc_cut = imc>.15;
%erode_imc= bwmorph(imc_cut, 'shrink', 1); figure(4); imagesc(erode_imc); colormap gray; axis square

%playing around with masking binary image and concolving
% figure(2)
% mask = imc_cut(95:120,10:20);
% mask = (mask+fliplr(mask)+flipud(mask)+fliplr(flipud(mask)))/4;
% 
% imagesc(mask); colormap gray;
% 
% imc_cutc = conv2(double(imc_cut)*256, mask, 'same'); 
% figure(3); imagesc(imc_cutc); colormap gray; axis square



%% find local maxima by making sure a pixel is greater than left, right, up, down
%% neighbors
% figure(2)
% imcD=imc([2:end 1],:);                                
% imcU=imc([end 1:end-1],:);                            
% imcR=imc(:,[2:end 1]);                     
% imcL=imc(:,[end 1:end-1]);                 
% imagesc(imc>=max(imcL,imcR));
% imagesc(imc>=max(imcU,imcD));
% imagesc(imc>=max(.1,max(max(imcL,imcR),max(imcU,imcD))));
%% unfortunately, a bunch of spurious points show up


%% square the correlation to make peaks more distinct
figure(3); imagesc(imc.^2 .* (imc>0)); axis square
imc2=imc.^2 .* (imc>0);  imc2=imc2/max(max(imc2));
imc2L=imc2([2:end 1],:);  
imc2R=imc2([end 1:end-1],:);                     
imc2U=imc2(:,[2:end 1]);    
imc2D=imc2(:,[end 1:end-1]);
figure(2); imagesc(imc2>=max(.1,max(max(imc2L,imc2R),max(imc2U,imc2D))));
%% much better!!!  (but not perfect)

%% let's make a more local mask... multiply by gaussian
maskM = fspecial('gaussian', size(mask2), 5);

%maskM=(ones(26,1)*exp(-(-14:14).^2/50)).*((exp(-(-27:27)'.^2/200))*ones(1,14));
figure(2); imagesc(maskM)
figure(2); imagesc(maskM.*mask2)
mask3=maskM.*mask2;
figure(3)
imc3 = conv2(im,mask3,'same'); imagesc(imc3); colormap gray; axis square
figure(3); imagesc(imc3.^2 .* (imc3>0)); axis square                           
figure(3); imagesc(imc .* imc3 .* (imc>0) .* (imc3>0)); axis square
imm = imc .* imc3 .* (imc>0) .* (imc3>0); imagesc(imm); axis square
imm=imm/(max(max(imm)));                                           
imagesc(imm>.065); axis square
%imagesc(imm); axis square
%% is this better?.... not clear....
% imm_cut = bwmorph(imm>.065, 'erode');
% figure(2); imagesc(imm_cut);colormap gray; axis square

%% make sure we can find the peaks of each blob... stick with imc
figure(1); imagesc(imc2); axis square
figure(2); imagesc(imc2>=max(imc2U,imc2D) & imc2>0); axis square
figure(3); imagesc(imc2>=max(imc2L,imc2R) & imc2>0); axis square; 
figure(2); imagesc(.5*(imc2>=max(imc2U,imc2D) & imc2>0)+imc2);
figure(3); imagesc(.5*(imc2>=max(imc2L,imc2R) & imc2>0)+imc2);
imc2c = (imc2>=max(max(imc2L,imc2R),max(imc2U,imc2D)) & imc2>.05);         %
figure(1); imagesc(imc2c + imc2);    pause;                             
%% yes, we get them... but sometimes we get a few extra peaks per blob

%% does it get any better if we first smooth the convolution image?
imc2a = conv2(imc2,ones(5)/25,'same');
imc2aL=imc2a([2:end 1],:);  
imc2aR=imc2a([end 1:end-1],:);                     
imc2aU=imc2a(:,[2:end 1]);    
imc2aD=imc2a(:,[end 1:end-1]);
imc2ac = (imc2a>=max(max(imc2aL,imc2aR),max(imc2aU,imc2aD)) & imc2a>.1);
figure(1); imagesc(imc2ac + imc2a)                                 
% no significant difference.  few points... could be good or bad

% find centers of DX... stick with imc2c
figure(1); imagesc(imc2c + imc2)                                 
[cI,cJ] = find(imc2c);
hold on; plot(cJ,cI,'go'); hold off
dI=13; dJ=7;
maskN=(ones(1+2*dI,1)*exp(-((-dJ:dJ).^2)/20)) .* ...
      ((exp(-((-dI:dI)'.^2)/80))*ones(1,1+2*dJ));
% replace each IJ with local maskN-weighted centroid of imc2
N=length(cI);
 for n=1:N
  if cI(n)-dI>=1 & cI(n)+dI<=512 & cJ(n)-dJ>=1 & cJ(n)+dJ<=512
   localI = (cI(n)-dI):(cI(n)+dI); localJ = (cJ(n)-dJ):(cJ(n)+dJ);
   local = imc2( localI , localJ ) .* maskN;
   cI(n) = sum(sum((localI' * ones(1,1+2*dJ)).*local))/sum(sum(local));
   cJ(n) = sum(sum((ones(1+2*dI,1)*localJ).*local))/sum(sum(local));
  end
 end
hold on; plot(cJ,cI,'r+'); hold off
%% new '+' locations are closer to blob centers than original 'o's

% then merge any points within (scaled) distance 10 of each other
merged = zeros(1,N);
 for n=1:N
  if ~merged(n)
   d = sqrt( (cI-cI(n)).^2/3.2 + (cJ-cJ(n)).^2 );
   nbr = find(d<10);
   if length(nbr)>1
     cI(n) = mean(cI(nbr));  cJ(n) = mean(cJ(nbr));
     merged(nbr(find(nbr~=n)))=1;
   end
  end
 end
cI=cI(find(merged==0)); cJ=cJ(find(merged==0));
hold on; plot(cJ,cI,'b*'); hold off
%% blue '*' are the final centers... most are pretty good

figure(2); imagesc(imc2); hold on; plot(cJ,cI,'r+'); hold off
figure(1); imagesc(im); hold on; plot(cJ,cI,'r+'); hold off
unhappy = ( cJ - cI/3 > 470 );
cI=cI(find(unhappy==0)); cJ=cJ(find(unhappy==0));
figure(2); imagesc(imc2); hold on; plot(cJ,cI,'r+'); hold off
figure(1); imagesc(im); hold on; plot(cJ,cI,'r+'); hold off


%%%% alternative to the above: find center-of-mass of each threshold blob
[cIp,cJp]=patchIJ(imc>.3*max(max(imc)));
imagesc(patchsize(imc>.3*max(max(imc)))); colormap hot
hold on; plot(cIp,cJp,'b.'); hold off
%% cI=cIp; cJ=cJp;  % uncomment to use these tile centers...
  
% now find four neighbors for each cell, if possible... check typical locations
figure(3); imagesc(mask2);
% expect neighbors at I+-23, J+-7
N=length(cI);
nUL=zeros(1,N); nUR=zeros(1,N); nLL=zeros(1,N); nLR=zeros(1,N); 
d=[0 0 0 0]; dN=0; dI=[0 0 0 0]; dJ=[0 0 0 0];
figure(2); imagesc(imc2); hold on; plot(cJ,cI,'r+'); 
 for n=1:N
   dUL = sqrt( (cI-(cI(n)-23)).^2/3.2 + (cJ-(cJ(n)-7)).^2 ); nbrUL = find(dUL<5);
   dUR = sqrt( (cI-(cI(n)-23)).^2/3.2 + (cJ-(cJ(n)+7)).^2 ); nbrUR = find(dUR<5);
   dLL = sqrt( (cI-(cI(n)+23)).^2/3.2 + (cJ-(cJ(n)-7)).^2 ); nbrLL = find(dLL<5);
   dLR = sqrt( (cI-(cI(n)+23)).^2/3.2 + (cJ-(cJ(n)+7)).^2 ); nbrLR = find(dLR<5);
   if length(nbrUL)>0, nUL(n) = nbrUL(1); else nUL(n)=0; end
   if length(nbrUR)>0, nUR(n) = nbrUR(1); else nUR(n)=0; end
   if length(nbrLL)>0, nLL(n) = nbrLL(1); else nLL(n)=0; end
   if length(nbrLR)>0, nLR(n) = nbrLR(1); else nLR(n)=0; end
   if nUL(n)*nUR(n)*nLL(n)*nLR(n)>0
     lI=cI( [n n n n; nUL(n) nUR(n) nLL(n) nLR(n)] );
     lJ=cJ( [n n n n; nUL(n) nUR(n) nLL(n) nLR(n)] );
     h=line(lJ, lI); set(h, 'Color', [1 1 0]*.7);
     dI = dI+ abs(diff(lI,1,1)); dJ = dJ + abs(diff(lJ,1,1));
     d = d+ sqrt( diff(lI,1,1).^2 /3.2 + diff(lJ,1,1).^2 ); dN=dN+1;
   end
 end
d/dN
dI/dN
dJ/dN
hold off
%% grid looks good superimposed on correlation image

%% draw grid on real image
figure(3); colormap gray; imagesc(im); hold on; plot(cJ,cI,'r+'); 
 for n=1:N
   if nUL(n)*nUR(n)*nLL(n)*nLR(n)>0
     lI=cI( [n n n n; nUL(n) nUR(n) nLL(n) nLR(n)] );
     lJ=cJ( [n n n n; nUL(n) nUR(n) nLL(n) nLR(n)] );
     h=line(lJ, lI); set(h, 'Color', [1 1 0]*.7);
   end
 end
hold off

pause;

%% get histograms of dI and dJ to neighbors
LRdata=[]; ULdata=[];
 for n=1:N
   if nUL(n)*nUR(n)*nLL(n)*nLR(n)>0
     dTOP= sqrt( (cI(nUL(n))-cI(nUR(n)))^2 + (cJ(nUL(n))-cJ(nUR(n)))^2 );
     dBOT= sqrt( (cI(nLL(n))-cI(nLR(n)))^2 + (cJ(nLL(n))-cJ(nLR(n)))^2 );
     LRdata=[LRdata dTOP dBOT];
     dLEFT = sqrt( (cI(nUL(n))-cI(nLL(n)))^2 + (cJ(nUL(n))-cJ(nLL(n)))^2 );
     dRIGHT= sqrt( (cI(nUR(n))-cI(nLR(n)))^2 + (cJ(nUR(n))-cJ(nLR(n)))^2 );
     ULdata=[ULdata dLEFT dRIGHT];
   end
 end
figure(2); hist(LRdata  /512 *250,20); xlabel('perpendicular to helix axis, nm')
figure(3); hist(ULdata  /512 *250,20); xlabel('parallel to helix axis, nm')


cd ~wrightc/