Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add weights to Patlak linear regression #924

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/include/stir/modelling/PatlakPlot.h
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,7 @@ class PatlakPlot : public RegisteredParsingObject<PatlakPlot, KineticModel>
float _time_shift; //!< Shifts the time to fit the timing of Plasma Data with the Projection Data.
bool _in_correct_scale; //!< Switch to scale or not the model_matrix to the correct scale, according to the appropriate scale factor.
bool _in_total_cnt; //!< Switch to choose the values of the model to be in total counts or in mean counts.
bool _assume_poisson_distribution; //!< Assume that image is Poisson distributed, and weight linear regression accordingly [1/(counts/int{Cp(t)}^2)]
std::string _blood_data_filename; //!< Name of file in which the input function is stored
PlasmaData _plasma_frame_data; //!< Stores the plasma data into frames for brain studies
std::string _time_frame_definition_filename; //!< name of file to get frame definitions
Expand Down
16 changes: 15 additions & 1 deletion src/modelling_buildblock/PatlakPlot.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ set_defaults()
this->_time_shift=0.;
this->_in_correct_scale=false;
this->_in_total_cnt=false;
this->_assume_poisson_distribution=false;
}

const char * const
Expand Down Expand Up @@ -191,7 +192,7 @@ PatlakPlot::apply_linear_regression(ParametricVoxelsOnCartesianGrid & par_image,
frame_num<=num_frames ; ++frame_num )
{
patlak_x[frame_num-1]=patlak_model_array[1][frame_num]/patlak_model_array[2][frame_num];
weights[frame_num-1]=1;
weights[frame_num-1]=1; // if (_assume_poisson_distribution) handled later in the main loop.
}
{ // Do linear_regression for each voxel // for k j i
float slope=0.F;
Expand Down Expand Up @@ -219,7 +220,19 @@ PatlakPlot::apply_linear_regression(ParametricVoxelsOnCartesianGrid & par_image,
// (remember, these are integrals over the time frame, not single values at discrete t)
for ( frame_num = starting_frame;
frame_num<=num_frames ; ++frame_num )
{
patlak_y[frame_num-1]=dyn_image[frame_num][k][j][i]/patlak_model_array[2][frame_num];
if(this->_assume_poisson_distribution)
{
weights[frame_num-1]=dyn_image[frame_num][k][j][i]/(patlak_model_array[2][frame_num]*patlak_model_array[2][frame_num]);
// if is zero, leave as is, zero weight?
if(weights[frame_num-1]){
weights[frame_num-1]=1/weights[frame_num-1];
}else{
weights[frame_num-1]=0.00001;
}
}
}
// Apply the regression to this pixel
linear_regression(y_intersection, slope,
chi_square,
Expand Down Expand Up @@ -335,6 +348,7 @@ initialise_keymap()
this->parser.add_key("Time Shift", &this->_time_shift);
this->parser.add_key("In total counts", &this->_in_total_cnt);
this->parser.add_key("In correct scale", &this->_in_correct_scale);
this->parser.add_key("Poisson distributed images", &this->_assume_poisson_distribution);
this->parser.add_key("Time Frame Definition Filename", &this->_time_frame_definition_filename);
this->parser.add_stop_key("end Patlak Plot Parameters");
}
Expand Down