How to store the corner points of harris in n*2 array

hello everyone, i use this algorithm :

I1=load('image.jpg')

 I =double(I1);
%****************************
axes(handles.axes2);
imshow(I1);
k = waitforbuttonpress; 
set(handles.text2, 'Visible', 'on');
point1 = get(gca,'CurrentPoint');  %button down detected
rectregion = rbbox;  %%%return figure units
point2 = get(gca,'CurrentPoint');%%%%button up detected
point1 = point1(1,1:2); %%%extract col/row min and maxs
point2 = point2(1,1:2);
lowerleft = min(point1, point2);
upperright = max(point1, point2); 
ymin = round(lowerleft(1)); %%%arrondissement aux nombrs les plus proches
ymax = round(upperright(1));
xmin = round(lowerleft(2));
xmax = round(upperright(2));
set(handles.text2, 'Visible', 'off');
%***********************************
Aj=6;
cmin=xmin-Aj; cmax=xmax+Aj; rmin=ymin-Aj; rmax=ymax+Aj;
min_N=20;max_N=100;
%%%%%%%%%%%%%%Intrest Points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma=2; Thrshold=20; r=6; disp=1;
dx = [-1 0 1; -1 0 1; -1 0 1]; % The Mask 
    dy = dx';
    %%%%%%
    Ix = conv2(I(cmin:cmax,rmin:rmax), dx, 'same');   
    Iy = conv2(I(cmin:cmax,rmin:rmax), dy, 'same');
    g = fspecial('gaussian',max(1,fix(6*sigma)), sigma); %%%%%%Gaussien Filter

    %%%%%
    Ix2 = conv2(Ix.^2, g, 'same');  
    Iy2 = conv2(Iy.^2, g, 'same');
    Ixy = conv2(Ix.*Iy, g,'same');
    %%%%%%%%%%%%%%
    k = 0.04;
    R11 = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;
    R11=(1000/max(max(R11)))*R11;
    R=R11;
    ma=max(max(R));
    sze = 2*r+1; 
    MX = ordfilt2(R,sze^2,ones(sze));
    R11 = (R==MX)&(R>Thrshold); 
    count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));

    loop=0;
    while (((count<min_N)|(count>max_N))&(loop<30))
        if count>max_N
            Thrshold=Thrshold*1.5;
        elseif count < min_N
            Thrshold=Thrshold*0.5;
        end

        R11 = (R==MX)&(R>Thrshold); 
        count=sum(sum(R11(5:size(R11,1)-5,5:size(R11,2)-5)));
        loop=loop+1;
    end

    R=R*0;
    R(5:size(R11,1)-5,5:size(R11,2)-5)=R11(5:size(R11,1)-5,5:size(R11,2)-5);
    [r1,c1] = find(R);
    PIP=[r1+cmin,c1+rmin];%%IP 

   %%%%%%%%%%%%%%%%%%%%Display

   Size_PI=size(PIP,1);
   for r=1: Size_PI
   I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)-2)=255;
   I(PIP(r,1)-2:PIP(r,1)+2,PIP(r,2)+2)=255;
   I(PIP(r,1)-2,PIP(r,2)-2:PIP(r,2)+2)=255;
   I(PIP(r,1)+2,PIP(r,2)-2:PIP(r,2)+2)=255;

   end
   axes(handles.axes2);
   imshow(uint8(I));

the algorithm detect the points and display them in black image, i want to get the positions of these points (x y) to use them in other algorithm (dbscan), because the later (dbscan) needs a array of 2 columns (one for x and other for y) of the points and the rows depending to the number of interest points extracted from harris algorithm.

 epsilon=3;
MinPts=10;
IDX=DBSCAN(X,epsilon,MinPts);
PlotClusterinResult(I2, IDX);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);

As you see X is the array. I tried the PIP variabl but the result was all of the points are in circle -->noise.