why cntrd function error

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

[Image Analyst]

The documentation doesn't seem to mention a cntrd() function. What toolbox is it in? Regardless I guess you're doing it incorrectly. What is "the piv experiment"? Where is that code obtained from?

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[승곤 유]

from youtube (63) Coding Workshop: Particle Image Velocimetry Part 1 - YouTube

[승곤 유]

The problem seems to be that the third value of the cnt should be output, but it is not

Answer 1

[the cyclist]

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.

Answer 2

[승곤 유]

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

Answer 3

[승곤 유]

Answer 4

[승곤 유]

Answer 5

[승곤 유]

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