why cntrd function error

8 views (last 30 days)
승곤 유
승곤 유 on 21 May 2022
Commented: 승곤 유 on 22 May 2022
Hello, I have a question while writing the cntrd function. If you use the function to analyze the piv experiment, the index value is supposed to come out, but it doesn't come out, so I'm asking why
  3 Comments
승곤 유
승곤 유 on 22 May 2022
The problem seems to be that the third value of the cnt should be output, but it is not

Sign in to comment.

Answers (5)

the cyclist
the cyclist on 21 May 2022
Edited: the cyclist on 21 May 2022
If you are talking about this particle tracking code, you might want to reach out to that team. You might not find very many people here using it.
But, it looks like the last update was in 2008, so you might not have much luck there.
You could try to post a small, specific section of the code that is not working and maybe folks here can help.

승곤 유
승곤 유 on 22 May 2022
I'm sorry. I didn't upload the code
close all
clear
clc
%이미지 불러오기
data_1=imread('synthetic00014.jpg'); %unit8형태
data_2=imread('synthetic00015.jpg');
data1=double(data_1); %double형태
data2=double(data_2);
%window size/roi
boxsize=16;
overlap=10;
spacing=boxsize-overlap;
Nrow=floor((size(data1,1)-boxsize)/spacing); %80
Ncol=floor((size(data1,2)-boxsize)/spacing); %97
%정수값을 얻기 위해 floor을 사용
%자기 상관분석을 위해 for loop를 사용하여 계산
%for문 돌릴때 roi를 지정해서 값을 제한해준다면 될꺼 같은데
for row=1:Nrow
for col=1:Ncol
%1. center of window
rowc=boxsize/2+spacing*(row-1);
colc=boxsize/2+spacing*(col-1);
%2. 하위 열/행
subrow=(rowc-boxsize/2+1):(rowc+boxsize/2);
subcol=(colc-boxsize/2+1):(colc+boxsize/2);
%위치마다의 데이터값을 정의
subdata1=data1(subrow,subcol);
subdata2=data2(subrow,subcol);
%정규화 : Z-점수 정규화 : (X-평균)/표준편차. 이상치를 잘 처리한다
subdata1=(subdata1-mean(subdata1(:)))/std(subdata1(:));
subdata2=(subdata2-mean(subdata2(:)))/std(subdata2(:));
%imshow(data1)
%rectangle('Position', [rowc colc 16 16],'EdgeColor','r')
%교차상관함수 계산
corr_func=fftshift(ifft2(fft2(subdata1).*conj(fft2(subdata2))))/size(subdata1,1)/size(subdata1,2);
corr_func=padarray(corr_func,[boxsize boxsize]);
%imagesc(corr_func)
pk=pkfnd(corr_func,0.1,boxsize/2);
%(a,최대상대강도,b) 최대상대강도는 프로그램 성능항상을 위해 바꿀수 있다
cnt=cntrd(corr_func,pk,5); %index값을 기반으로 max를 해야하나, 값 도출이 없다.
%5:pk 주변 5 픽셀내에서 검색
%(1,3,5,7~~)
%가장 밝은 피크의 위치를 찾고 해당 값을 저장
[value index]=max(cnt(:,1));
cnt_output=cnt(index,:);
%최종변수지정 : 중심에서의 변위값. cnt값 때문에 도출된 값이 이상함. cnt에서 밝기값이 index로 나와야 하나
%불가
dispx(row,col)=-1*cnt_output(:,1)-(boxsize/2+1+boxsize);
dispy(row,col)=-1*cnt_output(:,1)-(boxsize/2+1+boxsize);
end
end
%plot
x=(0:(size(dispx,2)-1))*spacing+boxsize/2;
y=(0:(size(dispx,1)-1))*spacing+boxsize/2;
[x2, y2]=meshgrid(x,y);
imagesc(data1);
colormap gray
hold on
quiver(x2,y2,dispx,dispy,20,'r')

승곤 유
승곤 유 on 22 May 2022

승곤 유
승곤 유 on 22 May 2022

승곤 유
승곤 유 on 22 May 2022
function out=pkfnd(im,th,sz)
% finds local maxima in an image to pixel level accuracy.
% this provides a rough guess of particle
% centers to be used by cntrd.m. Inspired by the lmx subroutine of Grier
% and Crocker's feature.pro
% INPUTS:
% im: image to process, particle should be bright spots on dark background with little noise
% ofen an bandpass filtered brightfield image (fbps.m, fflt.m or bpass.m) or a nice
% fluorescent image
% th: the minimum brightness of a pixel that might be local maxima.
% (NOTE: Make it big and the code runs faster
% but you might miss some particles. Make it small and you'll get
% everything and it'll be slow.)
% sz: if your data's noisy, (e.g. a single particle has multiple local
% maxima), then set this optional keyword to a value slightly larger than the diameter of your blob. if
% multiple peaks are found withing a radius of sz/2 then the code will keep
% only the brightest. Also gets rid of all peaks within sz of boundary
%OUTPUT: a N x 2 array containing, [row,column] coordinates of local maxima
% out(:,1) are the x-coordinates of the maxima
% out(:,2) are the y-coordinates of the maxima
%CREATED: Eric R. Dufresne, Yale University, Feb 4 2005
%MODIFIED: ERD, 5/2005, got rid of ind2rc.m to reduce overhead on tip by
% Dan Blair; added sz keyword
% ERD, 6/2005: modified to work with one and zero peaks, removed automatic
% normalization of image
% ERD, 6/2005: due to popular demand, altered output to give x and y
% instead of row and column
% ERD, 8/24/2005: pkfnd now exits politely if there's nothing above
% threshold instead of crashing rudely
% ERD, 6/14/2006: now exits politely if no maxima found
% ERD, 10/5/2006: fixed bug that threw away particles with maxima
% consisting of more than two adjacent points
%find all the pixels above threshold
%im=im./max(max(im));
ind=find(im > th);
[nr,nc]=size(im);
tst=zeros(nr,nc);
n=length(ind);
if n==0
out=[];
display('nothing above threshold');
return;
end
mx=[];
%convert index from find to row and column
rc=[mod(ind,nr),floor(ind/nr)+1];
for i=1:n
r=rc(i,1);c=rc(i,2);
%check each pixel above threshold to see if it's brighter than it's neighbors
% THERE'S GOT TO BE A FASTER WAY OF DOING THIS. I'M CHECKING SOME MULTIPLE TIMES,
% BUT THIS DOESN'T SEEM THAT SLOW COMPARED TO THE OTHER ROUTINES, ANYWAY.
if r>1 & r<nr & c>1 & c<nc
if im(r,c)>=im(r-1,c-1) & im(r,c)>=im(r,c-1) & im(r,c)>=im(r+1,c-1) & ...
im(r,c)>=im(r-1,c) & im(r,c)>=im(r+1,c) & ...
im(r,c)>=im(r-1,c+1) & im(r,c)>=im(r,c+1) & im(r,c)>=im(r+1,c+1)
mx=[mx,[r,c]'];
%tst(ind(i))=im(ind(i));
end
end
end
%out=tst;
mx=mx';
[npks,crap]=size(mx);
%if size is specified, then get ride of pks within size of boundary
if nargin==3 & npks>0
%throw out all pks within sz of boundary;
ind=find(mx(:,1)>sz & mx(:,1)<(nr-sz) & mx(:,2)>sz & mx(:,2)<(nc-sz));
mx=mx(ind,:);
end
%prevent from finding peaks within size of each other
[npks,crap]=size(mx);
if npks > 1
%CREATE AN IMAGE WITH ONLY PEAKS
nmx=npks;
tmp=0.*im;
for i=1:nmx
tmp(mx(i,1),mx(i,2))=im(mx(i,1),mx(i,2));
end
%LOOK IN NEIGHBORHOOD AROUND EACH PEAK, PICK THE BRIGHTEST
for i=1:nmx
roi=tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1))) ;
[mv,indi]=max(roi);
[mv,indj]=max(mv);
tmp( (mx(i,1)-floor(sz/2)):(mx(i,1)+(floor(sz/2)+1)),(mx(i,2)-floor(sz/2)):(mx(i,2)+(floor(sz/2)+1)))=0;
tmp(mx(i,1)-floor(sz/2)+indi(indj)-1,mx(i,2)-floor(sz/2)+indj-1)=mv;
end
ind=find(tmp>0);
mx=[mod(ind,nr),floor(ind/nr)+1];
end
if size(mx)==[0,0]
out=[];
else
out(:,2)=mx(:,1);
out(:,1)=mx(:,2);
end
end
function [out,out1,out2]=cntrd(im,mx,sz,interactive)
%cntrd: calculates the centroid of bright spots to sub-pixel accuracy.
% Inspired by Grier & Crocker's feature for IDL, but greatly simplified and optimized
% for matlab
% INPUTS:
% im: image to process, particle should be bright spots on dark background with little noise
% ofen an bandpass filtered brightfield image or a nice fluorescent image
%
% mx: locations of local maxima to pixel-level accuracy from pkfnd.m
%
% sz: diamter of the window over which to average to calculate the centroid.
% should be big enough
% to capture the whole particle but not so big that it captures others.
% if initial guess of center (from pkfnd) is far from the centroid, the
% window will need to be larger than the particle size. RECCOMMENDED
% size is the long lengthscale used in bpass plus 2.
%
%
% interactive: OPTIONAL INPUT set this variable to one and it will show you the image used to calculate
% each centroid, the pixel-level peak and the centroid
%
% NOTE:
% - if pkfnd, and cntrd return more then one location per particle then
% you should try to filter your input more carefully. If you still get
% more than one peak for particle, use the optional sz parameter in pkfnd
% - If you want sub-pixel accuracy, you need to have a lot of pixels in your window (sz>>1).
% To check for pixel bias, plot a histogram of the fractional parts of the resulting locations
% - It is HIGHLY recommended to run in interactive mode to adjust the parameters before you
% analyze a bunch of images.
%
%OUTPUT: a N x 3 array containing, x, y and brightness for each feature
% out(:,1) is the x-coordinates
% out(:,2) is the y-coordinates
% out(:,3) is the brightnesses
% out(:,4) is the sqare of the radius of gyration
%
%CREATED: Eric R. Dufresne, Yale University, Feb 4 2005
% 5/2005 inputs diamter instead of radius
% Modifications:
% D.B. (6/05) Added code from imdist/dist to make this stand alone.
% ERD (6/05) Increased frame of reject locations around edge to 1.5*sz
% ERD 6/2005 By popular demand, 1. altered input to be formatted in x,y
% space instead of row, column space 2. added forth column of output,
% rg^2
% ERD 8/05 Outputs had been shifted by [0.5,0.5] pixels. No more!
% ERD 8/24/05 Woops! That last one was a red herring. The real problem
% is the "ringing" from the output of bpass. I fixed bpass (see note),
% and no longer need this kludge. Also, made it quite nice if mx=[];
if nargin==3
interactive=0;
end
if sz/2 == floor(sz/2)
warning('sz must be odd, like bpass');
end
if isempty(mx)
warning('there were no positions inputted into cntrd. check your pkfnd theshold')
out=[];
return;
end
r=(sz+1)/2;
%create mask - window around trial location over which to calculate the centroid
m = 2*r;
x = 0:(m-1) ;
cent = (m-1)/2;
x2 = (x-cent).^2;
dst=zeros(m,m);
for i=1:m
dst(i,:)=sqrt((i-1-cent)^2+x2);
end
ind=find(dst < r);
msk=zeros([2*r,2*r]);
msk(ind)=1.0;
%msk=circshift(msk,[-r,-r]);
dst2=msk.*(dst.^2);
ndst2=sum(sum(dst2));
[nr,nc]=size(im);
%remove all potential locations within distance sz from edges of image
ind=find(mx(:,2) > 1.5*sz & mx(:,2) < nr-1.5*sz);
mx=mx(ind,:);
ind=find(mx(:,1) > 1.5*sz & mx(:,1) < nc-1.5*sz);
mx=mx(ind,:);
[nmx,crap] = size(mx);
%inside of the window, assign an x and y coordinate for each pixel
xl=zeros(2*r,2*r);
for i=1:2*r
xl(i,:)=(1:2*r);
end
yl=xl';
pts=[];
%loop through all of the candidate positions
for i=1:nmx
%create a small working array around each candidate location, and apply the window function
tmp=msk.*im((mx(i,2)-r+1:mx(i,2)+r),(mx(i,1)-r+1:mx(i,1)+r));
%calculate the total brightness
norm=sum(sum(tmp));
%calculate the weigthed average x location
xavg=sum(sum(tmp.*xl))./norm;
%calculate the weighted average y location
yavg=sum(sum(tmp.*yl))./norm;
%calculate the radius of gyration^2
rg=(sum(sum(tmp.*dst2))/ndst2);
%concatenate it up
%pts=[pts,[mx(i,1)+xavg-r,mx(i,2)+yavg-r,norm,rg]']; if u wanna find
%the total brightness and raius of gyration^2.
x_final = mx(i,1)+xavg-r ;
y_final = mx(i,2)+yavg-r ;
pts=[pts,[x_final,y_final]'];
%OPTIONAL plot things up if you're in interactive mode
if interactive==1
imagesc(tmp)
axis image
hold on;
plot(xavg,yavg,'x')
plot(xavg,yavg,'o')
plot(r,r,'.')
hold off
pause
end
end
out=pts';
end

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!