cabouman / mbircone

BSD 3-Clause "New" or "Revised" License
11 stars 9 forks source link

Add ROR control #89

Closed cabouman closed 2 years ago

cabouman commented 2 years ago

We need to add the following parameters to the interface of MBIRCONE: num_rows - controls number of rows in the reconstruction num_cols - controls number of columns in the reconstruciotn num_slices - controls number of slices in the reconstruction slice_offset - controls vertical offset of the center slice for the reconstruction

dyang37 commented 2 years ago

In the current MBIRCONE interface there's a variable called ror_radius, which specifies the horizontal radius of reconstruction in ALU. If both ror_radius and (num_rows, num_cols) are provided, do we want (num_rows, num_cols) overwrite ror_radius?

cabouman commented 2 years ago

Diyu, you are confused. Don't worry about ror_radius.

bommaritom commented 2 years ago

Hi all, is this correct:

num_rows - controls number of rows in the ROI num_cols - controls number of columns in the ROI num_slices - controls number of slices in the ROI slice_offset - controls vertical offset of the center slice for the ROI

If these control the ROR, then they will override the parameter ror_radius. This doesn't seem correct based on the discussion from earlier and I don't think this is what we want.

If they control the ROI (as described above) then there is no issue; these four will control the ROI and ror_radius controls the ROR. In that case, the ROR can only be controlled by modifying its radius (meaning it cannot be moved vertically) while the ROI can be freely moved based on the four parameters given.

dyang37 commented 2 years ago

Currently the implementation is assuming that (num_rows, num_cols, num_slices, slice_offset) is specifying the ROR (in branch ror_dev):

num_rows (int, optional): [Default=None] Integer number of rows in reconstructed image.
  If None, automatically set.
num_cols (int, optional): [Default=None] Integer number of columns in reconstructed image.
  If None, automatically set.
num_slices (int, optional): [Default=None] Integer number of slices in reconstructed image.
  If None, automatically set.
slice_offset(int, optional): [Default=0.0] slice offset in :math:`ALU` from the center of ROR. 

Given ROR defined by (num_rows, num_cols, num_slices, slice_offset), the new ROI should always be the intersection of old ROI and ROR.

cabouman commented 2 years ago

No, ror_radius is different. It effectively controls the radius of a binary mask that determines which pixels are updated. Pixels inside the ror_radius are updated. Pixels outside ror_radius are set to 0.

The variables num_rows, num_cols, num_slices, and slice_offset control the size of the array being reconstructed.

The two things are independent. The array size parameters are integers, but ror_radius is a float that has ALU units.

If the ror_radius is set to be larger than the maximum radius of the reconstruction array, then all the pixels will be reconstructed in the 3D array. However, this us usually not what you want.

bommaritom commented 2 years ago

In compute_img_size, the docstring says

- **ROR (list)**: Region of reconstruction that specifies the size of the reconstructed image. 
A list of 3 integer, [num_img_slices, num_img_rows, num_img_cols]. However, the valid region of interest (ROI) is a subset of ROR.

So, it looks like ROR controls a binary mask, and it also controls the actual image size outputted. We can also see this in MBIRModularUtilities.h:

struct ImageParams
{
    /* Location of the corner of the first voxel corresponding to
     (j_x, j_y, j_z) = (0, 0, 0). */
    float x_0;
    float y_0;
    float z_0;

    /* Number of voxels in x, y, z direction. */
    long int N_x;
    long int N_y;
    long int N_z;

    /* Dimensions of a voxel */
    float Delta_xy;
    float Delta_z;

    /**
     *      Region of Interest (roi) parameters
     */
    /* indices of the first voxels in the roi */
    long int j_xstart_roi;
    long int j_ystart_roi;
    long int j_zstart_roi;

    /* indices of the last voxesl in the roi */
    long int j_xstop_roi;
    long int j_ystop_roi;
    long int j_zstop_roi;

    long int N_x_roi;
    long int N_y_roi;
    long int N_z_roi;
};

If we then reference cone3D.py, we see that ror_radius is currently used to compute x_0, y_0, z_0, N_x, N_y, N_z, Delta_xy, Delta_z. What I'm saying is, it looks to me like ROR and "image size" are synonymous, while ROI is a subset of that image. In general (while I don't fully understand the C code), it looks like our flexibility is limited to the parameters in this struct.

Given this, I would assume that num_rows, num_slices, num_cols, slice_offset should be in correspondence with the parameters labeled "Region of Interest (roi) parameters" in the C code, while ror_radius corresponds to the first three groups.

bommaritom commented 2 years ago

Diyu also pointed out that j_xstart_roi, j_xstop_roi, N_x_roi, etc. are not actually used in any meaningful way by the C code. They are printed out to the terminal and are returned by compute_img_size function (for use in demo code, for e.g.) and they exist in several auxiliary C functions that are never called.

In fact, I added the following lines to recon directly after the call to compute_img_params, and the code ran with no issues:

imgparams['j_xstart_roi'] = -1
imgparams['j_ystart_roi'] = -1

imgparams['j_xstop_roi'] = -1
imgparams['j_ystop_roi'] = -1

imgparams['j_zstart_roi'] = -1
imgparams['j_zstop_roi'] = -1

imgparams['N_x_roi'] = -1
imgparams['N_y_roi'] = -1
imgparams['N_z_roi'] = -1

Conclusion: It does not look like the C code uses these parameters.

bommaritom commented 2 years ago

ROI_Unused_In_C.pptx I have cobbled together some slides in order to make sense of code structure and what these ROI parameters are actually doing. They don't seem to affect anything in the C code.

dyang37 commented 2 years ago

Thanks for the slides Marco! @Bouman, Charles @.***> Professor, should we discuss the ROR control issue during Lilly's meeting tomorrow? Or would you prefer to set up a separate meeting? I have some preliminary implementation in branch "ror_dev", but before that I think we should meet up and discuss what ROR control means. Regards, Diyu


From: Marco Bommarito @.> Sent: Wednesday, August 31, 2022 4:38 PM To: cabouman/mbircone @.> Cc: Diyu Yang @.>; Assign @.> Subject: Re: [cabouman/mbircone] Add ROR control (Issue #89)

---- External Email: Use caution with attachments, links, or sharing data ----

ROI_Unused_In_C.pptxhttps://github.com/cabouman/mbircone/files/9464839/ROI_Unused_In_C.pptx I have cobbled together some slides in order to make sense of code structure and what these ROI parameters are actually doing. They don't seem to affect anything in the C code.

— Reply to this email directly, view it on GitHubhttps://github.com/cabouman/mbircone/issues/89#issuecomment-1233397245, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AB2Q7OAHKNRO6J7ZNFGU5XTV367DRANCNFSM566D5JYQ. You are receiving this because you were assigned.Message ID: @.***>

cabouman commented 2 years ago

The interface should allow the user to set the following variables num_rows - integer that controls number of rows in the reconstruction num_cols - integer that controls number of columns in the reconstruction num_slices - integer that controls number of slices in the reconstruction slice_offset - float that controls vertical offset of the center slice for the reconstruction in units of ALU

If any of these are not set by the user, then they can be set in the following way: num_rows = np.round( num_det_columns*( (delta_pixel_detector/delta_pixel_image)/magnification ) ) num_cols = np.round( num_det_channels*( (delta_pixel_detector/delta_pixel_image)/magnification ) ) num_slices = np.round( num_det_rows*( (delta_pixel_detector/delta_pixel_image)/magnification ) ) slice_offset = 0.0

The variables should be set by a function named auto_image_size() as is done in the SVMBIR package.

After discussion, we determined that it is probably best to remove the argument ror_radius entirely from the python interface.

In addition, we need to change the function isInsideMask(long int i_1, long int i_2, long int N1, long int N2) so that it generates a circular mask rather than an elliptical mask. More specifically, it should do something like:

isInsideMask(long int i_1, long int i_2, long int N1, long int N2) N = max( N1, N2 ) if( i1^2 + i^2 <= N^2 ) Return(1) otherwise Return(0)

Functions that will need to be changed cone3D.recon() - The inputs and defaults need to be changed as described above cone3D.project() - For this function, the values of num_rows, num_cols, and num_slices are always available through the dimensions of the image. So in this case, the geometry parameters can be either given or default to the parameters associated size of the image. Notice that with the new defaults, the num_slices=num_det_rows, num_cols=num_det_channels. cone3D.compute_img_size() - I don't think we will need this function any more. cone3D.extract_roi_from_ror() - I don't think we will need this function any more.

old content Regarding the input ror_radius: First, we should probably change the name ror_radius to roi_radius. Do exactly what SVMBIR does. Allow the user to set roi_radius, but if it is left unset, then set the default value as:

roi_radius = delta_pixel_image * max( num_rows, num_cols)

gbuzzard commented 2 years ago

Perhaps we should include 'recon' as in: num_recon_rows - integer that controls number of rows in the reconstruction num_recon_cols - integer that controls number of columns in the reconstruction num_recon_slices - integer that controls number of slices in the reconstruction recon_slice_offset - float that controls vertical offset of the center slice for the reconstruction in units of ALU

or perhaps these could be encapsulated in a dict recon_shape['num_rows'] recon_shape['num_cols'] recon_shape['num_slices'] recon_shape['slice_offset']

cabouman commented 2 years ago

This is now complete. Hooray!