How to get a matrix which acquires a pixel values along a spiral trajectory of Image.

Hi, I'm working on compressive sensing for laser scanning applications. I want to sample the image along a continuous trajectory and get the pixel values. I'm able to get the pixel values along a continuous trajectory. In this case, a spiral trajectory. However, I don't know how to obtain the sampling matric C in y = Cx, where x is the input vector image, C is the sampling matrix which, on multiplying, gives values along the defined trajectory and y is the pixel values along that trajectory.

clear all;
clc;
%% Goal is to acuqire a pixel values along a spiral trajectory of Image (X).
% i.e y = Cx;
% where y is a m*1 measument matrix
% C is a m*N sampling matrix
% X is a N by 1 vectorized Image. 
% N = n1*n2 where n1, n2 are size of the image
% Paramters for Spiral
dSpirals = 30;
nSpirals = 20;
CompresRatio = 0.4; 
InputImage = imread('barbara1.png');
[height width] = size(InputImage);
N = height*width;
m = round(CompresRatio*N);
nPoints = m;
theta = linspace(0,360*nSpirals, nPoints);
% Define Spiral
x = round(((width/2)+1) +(theta/dSpirals).*cosd(theta));
y = round(((height/2)+1) +(theta/dSpirals).*sind(theta));
plot(x,y);
axis([0 width 0 height]);
grid on;
%% Measure the pixel values along the spiral trajectory
Measure = zeros(m,1);
for i = 1:m
    Measure(m) = InputImage(x(i),y(i));
end
%% How to find C in y = CX where:
% y = Measure
% X = InputImageV %i.e., InputImageV = InputImage(:);

[Matt J]

Are you sure you need the matrix? WHat's wrong with doing,

I=ind2sub(size(Image),x,y)

Measure=Image(I);

[Ajay Gunalan]

I need the matrix because in compressive sensing, we solve for x in y = Cx, and I need to pass both measure values y and the sampling matrix C to the reconstruction algorithm as seen in this line.

Compressive Sening Reconstruction GitHub Link

In this specific case, I just need the sampling matrix. But, for other reconstruction algorithms, I need:

measure values(y) = (Sampling Matrix C along a trajectory)*(transform domain)*(sparse x)

The goal here is to solve for x. So, Yes, I need the matrix to pass it to the reconstruction algorithm. Transform domain can be shearlet, wavelet, fourier, etc..

[Matt J]

I don't see how you'd be able to solve for x, because C is non-square. If C were square and invertible, it would be a permutation matrix. You don't need the matrix form of C to solve an equation which is just a permutation.

In any case, if you are certain that C is what you want, you will hopefully consider Accept-clicking my answer below.

[Ajay Gunalan]

Just as an optional comment, I'll add these. I was able to reconstruct the image (X) in spite of C being a non-square matrix. Click here for the result and code. The result is not perfect. However, I'm working on it to get the optimal result.

Answer 1

[Matt J]

sz=size(Image);
J=sub2ind(sz,x,y);
I=(1:numel(x))';

C=sparse(I,J,1);

[Ajay Gunalan]

I tried your answer as you can see from Matlab gitHub link for latest verison. (commit link) But, I get the following error when I check the SamplingMatrix C

Error using  * 
MTIMES (*) is not fully supported for integer classes. At least one argument must be scalar.

Error in spiral (line 41)
MeasureCheck = SamplerMat*InputImageV;

[Matt J]

You must conver the image to double floats.

[Ajay Gunalan]

Thanks the previous issue is solved. But, the sampler matrix C size doesn't match the size of InputImageV of size 262144. by 1. The sampler matrix has size 104858*249614. But it should be 104858*262144. Beacuse the image size is 512 (n1) by 512(n2). Thus X size is n1*n2 = 262144.

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches
the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.

Error in spiral (line 41)
MeasureCheck = SamplerMat*InputImageV; %y=Cx

[Matt J]

This should fix it.

sz=size(Image);
J=sub2ind(sz,x,y);
I=(1:numel(x))';

C=sparse(I,J,1,numel(x), prod(sz));

[Ajay Gunalan]

Thank you so much for solving it quickly. It solved my issue. I could have never solved this on my own. I still don't understand your solution. But it works. I will try to understand it. Once again, thanks a lot..!