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

cd ~winfree/Research/DNA/AFM_old/

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

pause;
    
    %%
    %% let's find a location that has a DX blob in the middle
    %% convolving with this (symmetrized) mask should locate the DX centers
    mask = im(263:343,152:232);
    mask = (mask+fliplr(mask)+flipud(mask)+fliplr(flipud(mask)))/4;     
    % assume mask/image is perfectly symmetric
    figure(2);
        imagesc(mask); 
        colormap gray; 
        axis square;
        title 'mask'

    imc = conv2(im,mask,'same'); 
    figure(3);
        imagesc(imc);  
        colormap gray; 
        axis square;
        title 'mask applied to raw image'
    %% convolution peaks aren't so sharp, but they are clear.
    %% mask is too large, encompassing several tiles 
    %%      this results in a convoluted image with false tiles created
    %%      where there were obvious blank spots in the original image.
    %% the larger mask also slows performance because of additional
    %% conputations required.
    %%

pause;

    %%
    %% but can we get a better mask from the convolution image?
    mask = imc(203:257,202:230);  
    mask = (mask+fliplr(mask)+flipud(mask)+fliplr(flipud(mask)))/4;
    % assume mask/image is perfectly symmetric
    figure(2);
        imagesc(mask);  
        colormap gray; 
        axis square;
        title 'mask from convulted image'

    imc = conv2(im,mask,'same');
    figure(3);
        imagesc(imc); 
        colormap gray; 
        axis square;
        title 'new convultion with second mask'
    %imagesc(imc>1000); 
    %% looks a bit better to me, maybe
    %% mask problem held over from previous section.
    %% this convultion and mask overwrite the previous ones.
    %%

pause;

    %%
    %% does it help to smooth it symmetrically?
    figure(2); 
        x = -10:10; 
        gs = (exp(-x.^2/4)/2-exp(-x.^2/16)/4); 
        plot(x,gs)   
        title 'laplacian filter'
        % matlab has programed laplacian filters.
    figure(3); 
        img = conv2(gs,gs,imc,'same'); 
        imagesc(img); 
        axis square;
        title 'convoluted image smoothed with gaussian filter'
    %% maybe looks worse.  Hmmm.
    %% mask problem heldover
    %% although this image shows the image defects very clearly (to the 
    %% human eye), such as row insertions
    %%

pause;

    %%
    %% or we could smooth the (already smooth) mask...
    figure(2)
        x = -10:10; 
        g = (exp(-x.^2/4)/2); 
        plot(x,g) 
        title 'gaussian filter'
        % matlab has a function for creating a gaussian filte rof arbitrary
        % size and standard deviation. 
    figure(3)
        mask2 = conv2(g,[1],mask,'same'); 
        imagesc(mask2); 
        axis square
        title 'gaussian smoothed mask'

pause;

    figure(2)
        imc = conv2(im,mask2,'same'); 
        imagesc(imc); 
        colormap gray; 
        axis square
        title 'convultion of raw image with gaussian smoothed mask'
    figure(3)
        imc = imc/max(max(imc));
        imagesc(im2bw(imc, graythresh(imc)));                      % THRESHOLDING performed here   
        title 'thresholded convoluted image'
    %% does look a tad better.
    %%
    
pause;

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

pause;
    
    %%
    %% square the correlation to make peaks more distinct
    figure(3); 
        imagesc(imc.^2 .* (imc>0)); 
        axis square
        % some of the fainter points are lost.
        title 'correlation squared'
    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))));
        title 'detected maxima'
    %% much better!!!  (but not perfect)
    %%

pause;

    %%
    %% let's make a more local mask... multiply by gaussian
    maskM=(ones(55,1)*exp(-(-14:14).^2/50)).*((exp(-(-27:27)'.^2/200))*ones(1,29));
    % there are much easier ways to make a gaussian...
    figure(2); 
        imagesc(maskM)
        title 'mask'
    figure(3); 
        mask3 = maskM .* mask2;
        imagesc(mask3)
        title 'gaussian filtered local mask'
    
pause;
    
    figure(2)
        imc3 = conv2(im,mask3,'same'); 
        imagesc(imc3); 
        colormap gray; 
        axis square
        title 'raw image convoluted with local mask'
    figure(3); 
        imagesc(imc3.^2 .* (imc3>0)); 
        axis square               
        colormap gray; 
        title 'convolution squared'
    figure(4); 
        imm = imc .* imc3 .* (imc>0) .* (imc3>0); 
        imagesc(imm); 
        axis square 
        colormap gray; 
        title 'previous and new convultion combined'
        
pause;
        
    figure(3)
        imm=imm/(max(max(imm)));                                           
        imagesc(imm>.1);                    % THRESHOLDING
        axis square
        colormap gray
        title 'thresholded convoltion image'
    %% is this better?.... not clear....
    %% yes, local features are not lost. crystal defects stay as defects
    %% rather than being falsely replaced by nice lattices.
    %%

pause;
    
    %%
    %% this appears to be repeated code. at the least repeated process on a
    %% different image.
    %% make sure we can find the peaks of each blob... stick with imc2 
    figure(1); 
        imagesc(imc2); 
        axis square
        title 'convoluted image'
    figure(2); 
        imagesc(imc2>=max(imc2U,imc2D) & imc2>0); 
        axis square
        title 'up down maxima'
    figure(3); 
        imagesc(imc2>=max(imc2L,imc2R) & imc2>0); 
        axis square
        title 'left right maxima'

pause;

    figure(2); 
        imagesc(.5*(imc2>=max(imc2U,imc2D) & imc2>0)+imc2);
        title 'up down combined with original'
    figure(3); 
        imagesc(.5*(imc2>=max(imc2L,imc2R) & imc2>0)+imc2);
        title 'left right combined with original'
    imc2c = (imc2>=max(max(imc2L,imc2R),max(imc2U,imc2D)) & imc2>.1);
    figure(1); 
        imagesc(imc2c + imc2)       
        title 'combine up down and left right'
    %% yes, we get them... but sometimes we get a few extra peaks per blob 
    %%

pause;

    %%
    %% 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)   
        title 'combined center detectors'
    %% no significant difference.  few points... could be good or bad
    %%

pause;

    %%
    %% find centers of DX... stick with imc2c
    figure(1); 
        imagesc(imc2c + imc2)  
        
    [cI,cJ] = find(imc2c);
    figure(2)
        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
    figure(3)
        hold on; 
        plot(cJ,cI,'r+'); 
        hold off
        
    %% new '+' locations are closer to blob centers than original 'o's
    %%

pause;

    %%
    %% 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));
    figure(2)
        hold on; 
        plot(cJ,cI,'b*'); 
        hold off
        
    %% blue '*' are the final centers... most are pretty good
    %%

pause;

    %% replot image and convoluted image with centers
    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
        title 'convoluted image with centers superimposed'
    figure(1); 
        imagesc(im); 
        hold on; 
        plot(cJ,cI,'r+'); 
        hold off
        title 'raw iamge with centers superimposed'
    %%

pause;
  
%%
%% now find four neighbors for each cell, if possible... check typical locations
figure(3); 
    imagesc(mask2);
    title 'mask'
% 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
    title 'convoluted image with centers and connections superimposed'
%% grid looks good superimposed on correlation image

%% draw grid on real image
figure(1); 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

%% 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')

figure(3); 


%% by gosh, this *is* faster, and equivalent
mcc -x avecell_fast
unitcell2 = avecell_fast(im,cI,cJ,nUL,nUR,nLL,nLR,35,20);

unitcellB = avecell_fast(im,cI,cJ,nUL,nUR,nLL,nLR,100,60);
imagesc(unitcellB); axis equal  
unitcellBS = ...
  unitcellB+fliplr(unitcellB)+flipud(unitcellB)+fliplr(flipud(unitcellB));
imagesc([unitcellB unitcellB]); axis equal
imagesc([unitcellBS unitcellBS]); axis equal


figure(4); imagesc(xcorr2(im)); axis square

pause;
close all;

cd ~wrightc