#include "StdAfx.h" #include "simultaneous_diffusion.h" #include <math.h> using namespace std; float randomNo(float k) {//to have sufficient precision rand() returns a short int (2B) x2 -> now returns a float (4B) return ( k * ((float)rand()+(float)rand()*RAND_MAX)/(RAND_MAX*(RAND_MAX+1)) ); //internal pseudo random gen is bad. replace by better one //return random_number_gen.randMT(k); } unsigned int randomNo_int() {//to have sufficient precision rand() returns a short int (2B) x2 -> now returns a float (4B) return rand()+rand()*RAND_MAX; //return random_number_gen.randMT(k); } double gamma(double x) {//returns an approx of the gamma function //http://www.rskey.org/gamma.htm double gamm; gamm = sqrt( x*sinh(1/x)+1/810/pow(x, 6) ); gamm *= x/exp((double)1); gamm = pow (gamm, x); gamm *= sqrt(6.283185307 / x); //2*Pi return gamm; } double gamma_pdf(double t, double alpha, double beta) {//probability density function of the gamma distribution //parameters k=gamma=alpha >0 and theta=beta >0 - expect alpha in ]2,3] (>2 => 2 inflexion points) and beta in [1,2] return ( (t/beta) * exp(-t/beta) / beta / gamma(alpha) ); } inline unsigned short GetSqrt_rw(unsigned short *sqrt_table_rw, unsigned int value) {// Uses precomputed table to get square root, unless entry is larger than table //unsigned short sqrt_table_rw [ sqrt_table_last_precomputed_rw+1 ]; if (value > sqrt_table_last_precomputed_rw) return ( (unsigned short)sqrt( (float)value ) ); else return ( sqrt_table_rw[value] ); }; //float mini(float input [], unsigned char l) //{//uses sort algorithm (multiple passes), probably slower than just looking for the minimum (1 pass) // //l size of input (# elements) // sort(input, input + l); //built-in sort algorithm // return input[0]; //} float mini(float input [], unsigned char l) {//should be faster than using sort algorithm (above) unsigned char index; float result; result = input[0]; //l size of input (# elements) for (index=1; index<l; index++) if (result > input[index]) result = input[index]; return result; } inline float coord_mod_w_float (float coord, unsigned short int w_coord) {//mirroring coordinates according to the cube w_dimensions int division; division = (int)abs(coord) / (int)w_coord; float remainder; //remainder = abs(coord) % float(w_coord); //remainder = abs(coord) - float(w_coord)*(float)division ; //fmod? remainder = fmod ( abs(coord), (float)w_coord); //printf ("fmod of 5.3 / 2 is %lf\n", fmod (5.3,2) ); //fmod of 5.3 / 2 is 1.300000 //fractpart = modf (pi , &intpart); //printf ("%lf = %lf + %lf \n", param, intpart, fractpart); //3.141593 = 3.000000 + 0.141593 //if (division % 2) //division is odd // return (float)w_coord - remainder; //else //division is even // return remainder; if (division % 2) //division is odd remainder = (float)w_coord - remainder; if (remainder == (float)w_coord) remainder -= (float).0001; //ok up to dimensions < 1000 return remainder; } inline float coord_non_mod_float(float coord, unsigned short int w_coord) {//NOT mirroring coordinates UNREFERENCED_PARAMETER(w_coord); //to get rid of compiler warning message return coord; } inline unsigned short int coord_mod_w_int (int coord, unsigned short int w_coord) {//mirroring coordinates according to the cube w_dimensions int division; division = abs(coord) / (int)w_coord; unsigned short int remainder; //remainder = abs(coord) % float(w_coord); //remainder = abs(coord) - float(w_coord)*(float)division ; //fmod? remainder = (unsigned short int) (abs(coord) % w_coord); if (division % 2) //division is odd { remainder = w_coord - remainder;// - 1; if (coord > 0) //and positive coord -> need to shift index for proper voxel location remainder -- ; else if (remainder == w_coord) remainder -- ; } else //division is even if ((remainder != 0) && (coord < 0)) //and negative coord -> need to shift index for proper voxel location remainder -- ; return remainder; } inline unsigned short int coord_non_mod_int (int coord, unsigned short int w_coord) {//NOT mirroring coordinates UNREFERENCED_PARAMETER(w_coord); //to get rid of compiler warning message return (unsigned short int)coord; } //float maxi(float input [], unsigned char l) //{ // //l size of input (# elements) // sort(input, input + l); //built-in sort algorithm // return input[l-1]; //} sg_status_rw ValidateOuterSphereLayer_rw(DS_1b_cube* input_cube, unsigned short *sqrt_table_rw, unsigned short z, unsigned short y, unsigned short x, radius_map_type_rw prev_rad, radius_map_type_rw rad) {//from Vasa's code - sphere growing //function for finding a radius at a given spot... uses neighbors already calc radius values to speed things up short int dz, dy; int z_m_dz, z_p_dz, y_m_dy, y_p_dy, x_m_dx, x_p_dx; short int dx, dx_inner_lim, dx_outer_lim; unsigned int dz2, dz2_dy2, rad2 = rad*rad, rad_prev2 = prev_rad*prev_rad; for (dz = 0; dz <= (short int)rad; dz++) { dz2 = dz*dz; z_p_dz = z+dz; z_m_dz = z-dz; for (dy = 0; dy <= GetSqrt_rw(sqrt_table_rw, rad2 - dz2); dy++) { //dy^2 dz2_dy2 = dz2 + dy*dy; y_p_dy = y+dy; y_m_dy = y-dy; dx_outer_lim = GetSqrt_rw(sqrt_table_rw, rad2 - dz2_dy2); dx_inner_lim = ( dz2_dy2 > rad_prev2 ) ? (-1) : ( GetSqrt_rw(sqrt_table_rw, rad_prev2 - dz2_dy2) ); for (dx = dx_outer_lim; dx > dx_inner_lim; dx--) { x_p_dx = x+dx; x_m_dx = x-dx; //pos z, pos y if (z_p_dz < input_cube->wZ) { if (y_p_dy < input_cube->wY) { if ( (x_p_dx < input_cube->wX) && (input_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx)) ) return SG_FIBER_HIT_rw; if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && input_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx)) return SG_FIBER_HIT_rw; } if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) ) {//neg y if ( (x_p_dx < input_cube->wX) && (input_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx)) ) return SG_FIBER_HIT_rw; if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && input_cube->get_spot((unsigned short)z_p_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx)) return SG_FIBER_HIT_rw; } } //neg z if ( (z_m_dz != z_p_dz) && (z_m_dz >= 0) ) { if (y_p_dy < input_cube->wY) { if ( (x_p_dx < input_cube->wX) && (input_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_p_dx)) ) return SG_FIBER_HIT_rw; if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && input_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_p_dy, (unsigned short)x_m_dx)) return SG_FIBER_HIT_rw; } if ( (y_m_dy != y_p_dy) && (y_m_dy >= 0) ) {//neg y if ( (x_p_dx < input_cube->wX) && (input_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_p_dx)) ) return SG_FIBER_HIT_rw; if ( (x_m_dx != x_p_dx) && (x_m_dx >= 0) && input_cube->get_spot((unsigned short)z_m_dz, (unsigned short)y_m_dy, (unsigned short)x_m_dx)) return SG_FIBER_HIT_rw; } } } } } return SG_NO_ERROR_rw; } radius_map_type_rw GetExactRadius_rw(DS_1b_cube* input_cube, DSt_cube<radius_map_type_rw> *distMap, unsigned short *sqrt_table_rw, unsigned short z, unsigned short y, unsigned short x) {//from Vasa's code - sphere growing //function for checking if there are any intersections with fibers and outside sphere wall radius_map_type_rw prev_max_rad = 0; // Find the maximum neighboring radius from previous locations // + check not a fiber voxel (one single radius map cube for both!!) //done only for spot (xyz)-1, because of the scan order from the loop in the main() if ( (x > 0) && (prev_max_rad < distMap->datac(z, y, x-1)) && (!input_cube->get_spot(z, y, x-1)) ) prev_max_rad = distMap->datac(z, y, x-1); if ( (y > 0) && (prev_max_rad < distMap->datac(z, y-1, x)) && (!input_cube->get_spot(z, y-1, x)) ) prev_max_rad = distMap->datac(z, y-1, x); if ( (z > 0) && (prev_max_rad < distMap->datac(z-1, y, x)) && (!input_cube->get_spot(z-1, y, x)) ) prev_max_rad = distMap->datac(z-1, y, x); // If the three previous radius are all '0', then we are either on a boundary or one spot away... max radius is 1 // First try to verify the same size as greatest last space // -if last was '0', then no point // If previous max radius was 0, find maximum we can fit if ( prev_max_rad == 0 ) { while ( (ValidateOuterSphereLayer_rw(input_cube, sqrt_table_rw, z, y, x, prev_max_rad, prev_max_rad+1) != SG_FIBER_HIT_rw) && (prev_max_rad<254)) //<254 added in the case of an 8-bit resolution distMap -> avoid overflow prev_max_rad++; return (prev_max_rad); } // -if we fail the validation, radius is 1 smaller than last time else if (ValidateOuterSphereLayer_rw(input_cube, sqrt_table_rw, z, y, x, prev_max_rad-1, prev_max_rad) == SG_FIBER_HIT_rw) return (prev_max_rad - 1); //first make sure we don't create a sphere bigger than max allowed if ( prev_max_rad >= 254 ) //8bits case { char txt[1024]; sprintf(txt, "We might have a sphere bigger than max allowed radius (254) at (x,y,z)[%dx, %dy, %dz]", x, y, z); write_to_log(txt); return 254; //prev_max_rad; } // Next we try one bigger than last radius // -if we succeed, radius is one bigger than last // -if we fail, radius is same as last if ( ValidateOuterSphereLayer_rw(input_cube, sqrt_table_rw, z, y, x, prev_max_rad, prev_max_rad+1) == SG_FIBER_HIT_rw ) return (prev_max_rad); return (prev_max_rad + 1); } void simultaneous_diffusion(DS_1b_cube* input_cube, void * dlg_share) { //variables DSt_cube<radius_map_type_rw> *distMap = new DSt_cube<radius_map_type_rw>(input_cube->wZ, input_cube->wY, input_cube->wX); LinkedList<distance_structure> *distance_structure_list_ptr = new LinkedList<distance_structure> (1000000); //~23.1MB mem chunk //get all the parameter values from the dialog RW_Dialog* dlg = (RW_Dialog*) dlg_share; //initialize file names and directories char input_directory[1024], input_full_name[1024], work_directory[1024], output_full_name[1024]; { if ( 1 == dlg->CHECK_RW_BATCH_MODE.GetCheck() ) {//batch mode case,need to update "input_file_name_cut" dlg->input_file_name_cut = dlg->input_file_name; int index; index = dlg->input_file_name_cut.GetLength() - 1 ; while (index >= 0 ) { if ( dlg->input_file_name_cut[index] == '(' ) break; index--; } //if (left_bracket_indx > 0) // input_file_name_cut.Delete(left_bracket_indx, input_file_name.GetLength() ); if (index > 0) dlg->input_file_name_cut.Delete(index, dlg->input_file_name_cut.GetLength() ); } sprintf(input_directory, "%s", dlg->input_directory_name); sprintf(input_full_name, "%s%s", dlg->input_directory_name, dlg->input_file_name_cut); //if ( 1 == dlg->CHECK_RW_BATCH_MODE.GetCheck() ) // sprintf(work_directory, "%s%s", dlg->field_directory_work_str, "batch_mode\\"); //need to check if folder exists or create it //else sprintf(work_directory, "%s", dlg->field_directory_work_str); if ( 1 == dlg->CHECK_RW_BATCH_MODE.GetCheck() ) { //sprintf(output_full_name, "%s%s", work_directory, dlg->input_file_name_cut); char temp_field_name_output_str[1024]; sprintf(temp_field_name_output_str, "%s%s", dlg->input_file_name_cut, "_batch_mode"); dlg->field_name_output_str = temp_field_name_output_str; //if batch mode then output file name has to reflect input name //dlg->field_name_output_str = dlg->input_file_name_cut; } sprintf(output_full_name, "%sRandomWalk\\%s", dlg->field_directory_work_str, dlg->field_name_output_str); //write_to_log("Input file selected: %s", input_full_name); write_to_log("Output: %s", output_full_name); } //initialise local variables from dialog values bool mirror_coordinates, infinite_mirroring, save_full_RW_data, random_starting_point, sub_ROI_is_sphere, check_cv; float max_displacement_diag_pixels_sq; { mirror_coordinates = ( 1 == dlg->CHECK_RW_MIRROR.GetCheck() ); infinite_mirroring = ( 1 == dlg->CHECK_RW_MIRROR_INFINITE.GetCheck() ); save_full_RW_data = ( 1 == dlg->CHECK_SAVE_RW_RAW_DATA.GetCheck() ); random_starting_point = ( 1 == dlg->CHECK_RW_RANDOM_START.GetCheck() ); sub_ROI_is_sphere = ( 1 == dlg->CHECK_SUB_ROI_SPHERE.GetCheck() ); sub_ROI_is_sphere |= mirror_coordinates; //if mirror then force selection "sphere type of sub ROI" check_cv = ( 1 == dlg->CHECK_CV.GetCheck() ); max_displacement_diag_pixels_sq = pow((float)dlg->field_rw_MAX_DISPLACEMENT_DIAG_COEFF_float, 2) * \ ( pow((float)input_cube->wZ, 2) + pow((float)input_cube->wY, 2) + pow((float)input_cube->wX, 2) ); //in pixels } //initialise run variables bool max_distance_exceeded, rw_converged, rod_zero_aborted_enough_data_for_stats; unsigned int max_distance_exceeded_count; { max_distance_exceeded = false; max_distance_exceeded_count = 0; rw_converged = false; rod_zero_aborted_enough_data_for_stats = true; } //pointers to coord translation function - mirroring case float (*coord_trans_float) (float coord, unsigned short int w_coord); unsigned short int (*coord_trans_int) (int coord, unsigned short int w_coord); switch (mirror_coordinates) { case false: coord_trans_float = coord_non_mod_float; coord_trans_int = coord_non_mod_int; break; //case true: default: //to avoid Warning C4701: potentially uninitialized local variable 'coord_trans_int' used compiler warning message coord_trans_float = coord_mod_w_float; coord_trans_int = coord_mod_w_int; break; } float PIXEL_SIZE, PIXEL_SIZE_sq; unsigned int total_particle_per_roi; //# of tracers in each ROI unsigned short int total_file_number; //# of ROI(sub-cube) to select in the volume unsigned long int total_particle_end_of_run = 0; //# of tracers at the end of the simulation (for log) float dataCollectionPoint_increment; //in MICRONS!! float max_allowed_distance_pix; float rod, inv_sqrt_rod, sqrt_rod; //ratio of diffusivity float lambda_pix, lambdaForFiber_pix, lambda_micron, lambdaForFiber_micron; float max_jump_size; //number of times lambda_pix bar to limit jump size (set it to large value -> no limit) float max_jump_size_lambda_pix; //in pixels float PIXEL_SIZE_sq_over_2lambda_micron; float PIXEL_SIZE_sq_over_2lambdaForFiber_micron; float FPT_boundary_layer_coef; float FPT_boundary_layer_lambda_pix, FPT_boundary_layer_lambda_pix_4_fifth; //in pixels float FPT_boundary_layer_lambdaForFiber_pix, FPT_boundary_layer_lambdaForFiber_pix_4_fifth=0; //in pixels float two_FPT_boundary_layer_lambda_pix, two_FPT_boundary_layer_lambdaForFiber_pix; //in pixels float sq_MAX_ROI_RADIUS; float distance; float max_distance, max_displ2, max_displ2x, max_displ2y, max_displ2z; float D_D0x, D_D0y, D_D0z; float D_D0x_previous, D_D0y_previous, D_D0z_previous; unsigned long cv_num_tracers; float cv_criteria;//, cv_criteria_sq; float cv_percent_new_data; unsigned long long count_distances; unsigned long long num_jumps_pore_fpt, num_jumps_pore_nonfpt, num_jumps_fiber_fpt, num_jumps_fiber_nonfpt; float total_ROI_porosity; float coef_rod_x, coef_rod_y, coef_rod_z, coef_rod_sex, coef_rod_sey, coef_rod_sez, coef_rod0_x, coef_rod0_y, coef_rod0_z; float xROImin, xROImax, yROImin, yROImax, zROImin, zROImax; //ROI (sub-cube) coordinates int num_coeff; unsigned short int file_number; unsigned int min_required_num_per_bucket; //dependent on max pop bucket, see below int total_count_tracer_start_hit_pore, total_count_tracer_start_hit; { PIXEL_SIZE = dlg->field_pixel_size_float; //(float)1; PIXEL_SIZE_sq = pow(PIXEL_SIZE, 2); total_particle_per_roi = dlg->field_rw_MAX_NUM_TRACERS_int; //50000; total_file_number = (unsigned short int)dlg->field_rw_NUM_SUB_CUBES_int; //1; dataCollectionPoint_increment = dlg->field_rw_DATACOLLECTIONPOINT_INCREMENT_float; //10; max_allowed_distance_pix = dlg->field_rw_MAX_DISTANCE_float/(float)PIXEL_SIZE;//(float)1000000/(float)PIXEL_SIZE; //(1m -> pixels) //( rod=0 means pore diffusion only (no diff in fiber)) //-> need to take into account porosities of sub-cube //( rod=1 means open space diffusion rod = dlg->field_rw_ROD_float; //(float)0.04; //ratio of diffusivity //(usually between 0.02 to 0.8 which corresponds to RH 0% to 100% if (rod == 0) { sqrt_rod = -1; //this way tracer will NOT enter fiber space inv_sqrt_rod = 0; } else { sqrt_rod = sqrt(rod); //not to calculate it everytime inv_sqrt_rod = 1/sqrt_rod; } lambda_micron = dlg->field_rw_LAMBDA_float; lambda_pix = lambda_micron/(float)PIXEL_SIZE; //convert to pixels lambdaForFiber_micron = dlg->field_rw_LAMBDA_FOR_FIBER_float; lambdaForFiber_pix = lambdaForFiber_micron/(float)PIXEL_SIZE; //convert to pixels //float lambda_sq, lambdaForFiber_sq; //not to calculate them everytime //lambda_sq = pow(lambda, 2); //lambdaForFiber_sq = pow(lambdaForFiber, 2); max_jump_size = dlg->field_rw_MAX_JUMP_SIZE_float;//(float)500000; max_jump_size_lambda_pix = max_jump_size * lambda_pix; PIXEL_SIZE_sq_over_2lambda_micron = PIXEL_SIZE_sq / (2*lambda_micron); PIXEL_SIZE_sq_over_2lambdaForFiber_micron = PIXEL_SIZE_sq / (2*lambdaForFiber_micron); FPT_boundary_layer_coef = dlg->field_rw_FPT_BOUNDARY_LAYER_LAMBDA_float; FPT_boundary_layer_lambda_pix = FPT_boundary_layer_coef*lambda_pix; FPT_boundary_layer_lambda_pix_4_fifth = (float)4/(float)5*FPT_boundary_layer_lambda_pix; two_FPT_boundary_layer_lambda_pix = (float)2*FPT_boundary_layer_lambda_pix ; FPT_boundary_layer_lambdaForFiber_pix = FPT_boundary_layer_coef*lambdaForFiber_pix; FPT_boundary_layer_lambdaForFiber_pix_4_fifth = (float)4/(float)5*FPT_boundary_layer_lambdaForFiber_pix_4_fifth; two_FPT_boundary_layer_lambdaForFiber_pix = (float)2*FPT_boundary_layer_lambdaForFiber_pix; sq_MAX_ROI_RADIUS = pow((float)dlg->field_rw_SIZE_SUB_CUBES_float/2, 2); //in pixels if (infinite_mirroring) //sq_MAX_ROI_RADIUS = min(sq_MAX_ROI_RADIUS, max_displacement_diag_pixels_sq); //take the value from the dlg and calculated above for the inf mirroring case sq_MAX_ROI_RADIUS = max_displacement_diag_pixels_sq; max_distance = max_displ2 = max_displ2x = max_displ2y = max_displ2z = (float)0; D_D0x_previous = 0; D_D0y_previous = 0; D_D0z_previous = 0; //if (check_cv) //commented out to avoid warning not initialised error message { cv_num_tracers = dlg->field_rw_CV_NUM_TRACERS_int; cv_criteria = dlg->field_rw_CV_CRITERIA_float; //cv_criteria_sq = pow(cv_criteria, 2); cv_percent_new_data = dlg->field_rw_CV_NEW_DATA_float; } count_distances = 0; total_ROI_porosity = (float)0; xROImin = 0; xROImax = 0; yROImin = 0; yROImax = 0; zROImin = 0; zROImax = 0; total_count_tracer_start_hit_pore = 0; //used to get porosity from where (coord) the tracers land at their creation total_count_tracer_start_hit = 0; } //display message in RW dialog window char win_txt [1024]; sprintf(win_txt, "Simultaneous Diffusion: rod= %f ", rod); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, 2 * input_cube->wZ); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //file to save future results into (text) char input_base_name[1024]; //directories that are needed sprintf(input_base_name, "%s%s", input_directory, dlg->input_file_name_cut);//dlg->field_name_output_str); //init file base name char init_file_name[1024]; sprintf(init_file_name, "%sRandomWalk\\%s", work_directory, dlg->input_file_name_cut);//dlg->field_name_output_str); char init_file_out_name[1024]; sprintf(init_file_out_name, "%sRandomWalk\\%s", work_directory, dlg->field_name_output_str); char pre_file_out [1024]; sprintf(pre_file_out, "%s(%d_%d_%d).stack.cbin", init_file_out_name, input_cube->wX, input_cube->wY, input_cube->wZ); char init_file_name_for_radius[1024]; sprintf(init_file_name_for_radius, "%sRandomWalk\\%s", work_directory, dlg->input_file_name_cut); char radius_file_out [1024]; sprintf(radius_file_out, "%s_RADIUS(%d_%d_%d).u8.raw", init_file_name_for_radius, input_cube->wX, input_cube->wY, input_cube->wZ); //file for saving full RW raw data ofstream rw_raw_data_fout; //if (save_full_RW_data) //{ // char rw_raw_data_filename[1024]; // char* rw_raw_data_file; // sprintf(rw_raw_data_filename, "%s_rw_raw_data.csv", output_full_name);//init_file_name_for_radius); // rw_raw_data_file = rw_raw_data_filename; // // rw_raw_data_fout.open(rw_raw_data_file, ios::out | ios::app); // if( !rw_raw_data_fout.is_open() ) // { // char text[1024]; // sprintf(text,"Simultaneous Diffusion: error opening file for saving rw_raw_data: %s", rw_raw_data_filename); // write_to_log(text); // return; // } // //rw_raw_data_fout << "test" << endl; //} //file for saving CV raw data ofstream rw_cv_data_fout; ofstream rw_cv_data_detailed_fout; if (check_cv) { //rw_cv_data_fout char rw_cv_data_filename[1024]; char* rw_cv_data_file; sprintf(rw_cv_data_filename, "%s_rw_cv_data.csv", output_full_name);//init_file_name_for_radius); rw_cv_data_file = rw_cv_data_filename; rw_cv_data_fout.open(rw_cv_data_file, ios::out | ios::app); if( !rw_cv_data_fout.is_open() ) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving rw_cv_data: %s", rw_cv_data_filename); write_to_log(text); return; } //write header rw_cv_data_fout << "Input file selected, " << input_full_name <<endl; rw_cv_data_fout << "Output, " << output_full_name <<endl; rw_cv_data_fout << "saving rw_cv_data, " << rw_cv_data_filename <<endl; rw_cv_data_fout << "#tracers, D_D0x, D_D0y, D_D0z, %x, %y, %z, R2_x, R2_y, R2_z, cv_criteria <," << cv_criteria << endl; //rw_cv_data_detailed_fout char rw_cv_data_detailed_filename[1024]; char* rw_cv_data_detailed_file; sprintf(rw_cv_data_detailed_filename, "%s_rw_cv_data_detailed.csv", output_full_name);//init_file_name_for_radius); rw_cv_data_detailed_file = rw_cv_data_detailed_filename; rw_cv_data_detailed_fout.open(rw_cv_data_detailed_file, ios::out | ios::app); if( !rw_cv_data_detailed_fout.is_open() ) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving rw_cv_data: %s", rw_cv_data_detailed_filename); write_to_log(text); return; } //write header rw_cv_data_detailed_fout << "Input file selected, " << input_full_name <<endl; rw_cv_data_detailed_fout << "Output, " << output_full_name <<endl; rw_cv_data_detailed_fout << "saving rw_cv_data, " << rw_cv_data_detailed_filename <<endl; rw_cv_data_detailed_fout << "#tracers, D_D0x, D_D0y, D_D0z, %x, %y, %z, R2_x, R2_y, R2_z, cv_criteria <," << cv_criteria << endl; rw_cv_data_detailed_fout << ", av_dist, av_displ2X, av_displ2Y, av_displ2Z, count, "; rw_cv_data_detailed_fout << endl << endl; } //calculate distance maps, if file does not exist already if ( distMap->readFromFile(radius_file_out, (unsigned short)0) == DST_FILE_IO_ERROR) { write_to_log("Simultaneous Diffusion - creating distMap"); //precompute all possible sqrts (cuts run time in half), they are used many many times sprintf(win_txt, "Precomputing square root table..."); dlg->txt_field_progress.SetWindowText(win_txt); unsigned short sqrt_table_rw [ sqrt_table_last_precomputed_rw+1 ]; register unsigned int i; for (i = 0; i <= sqrt_table_last_precomputed_rw; i++) sqrt_table_rw[i] = (unsigned short) sqrt( (float)i ); unsigned short int x, y, z; sprintf(win_txt, "Distance Map 1/2 (pore space)"); dlg->txt_field_progress.SetWindowText(win_txt); // Clear distMap cube distMap->clear(); //do pore space map//DistanceMapPore (); for (z = 0; z < input_cube->wZ; z++) {// Get exact radius at each pore voxel, uses neighbor info to speed things up for (y = 0; y < input_cube->wY; y++) for (x = 0; x < input_cube->wX; x++) if ( !input_cube->get_spot(z, y, x) ) distMap->datac(z, y, x, GetExactRadius_rw(input_cube, distMap, sqrt_table_rw, z, y, x) ); //update progress bar dlg->pbar_simulation_progress.StepIt(); } write_to_log("distMap 1/2 (pore space) - done"); //if user requested termination of simulation if (dlg->stop_request) goto skip_fiber_space_distMap; //invert cube input_cube->invert(); sprintf(win_txt, "Distance Map 2/2 (fiber space)"); dlg->txt_field_progress.SetWindowText(win_txt); //do pore fiber map//DistanceMapFiber(); for (z = 0; z < input_cube->wZ; z++) {// Get exact radius at each fiber voxel, uses neighbor info to speed things up for (y = 0; y < input_cube->wY; y++) for (x = 0; x < input_cube->wX; x++) if ( !input_cube->get_spot(z, y, x) ) distMap->datac(z, y, x, GetExactRadius_rw(input_cube, distMap, sqrt_table_rw, z, y, x)); //update progress bar dlg->pbar_simulation_progress.StepIt(); } write_to_log("distMap 2/2 (fiber space) - done"); //invert cube back to original input_cube->invert(); //delete [] sqrt_table_rw; //save the cube file sprintf(radius_file_out, "%s_RADIUS", init_file_name); distMap->writeToFile(radius_file_out); write_to_log("distMap saved to file: %s", init_file_name_for_radius); } else write_to_log("Simultaneous Diffusion - already existing distMap loaded"); skip_fiber_space_distMap: char stats_results_file_out[1024]; char* stats_results_file; //stats results file sprintf(stats_results_file_out, "%s.simult_diff_results.csv", pre_file_out); stats_results_file = stats_results_file_out; ofstream stats_results_fout; stats_results_fout.open(stats_results_file, ios::out | ios::app);// | std::ios::binary | ios_base::right); if(!stats_results_fout.is_open()) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving stats results: %s", stats_results_file_out); write_to_log(text); return; } MTRand random_number_gen; //automatically seeds the MersienneTwister random gen //double j = random_number_gen.randMT(1000000.0); //{//PRNGs testing for evaluation via diehard test - g marsaglia // //seed the pseudorandom-number generator before calling rand. // srand( ( (unsigned)(time(NULL)*(double)(input_cube->wZ%211)*(double)(input_cube->wY%277)/(double)(input_cube->wX%293)) )% RAND_MAX ); //211, 277 & 293 prime // //for debug, test RandomNum generators // //for (file_number=1 ; file_number <= 65000; file_number++) // total_particle_per_roi = (unsigned int)( 3*pow(10.0,5.0) ); // write_to_log("start %d randomNo", total_particle_per_roi); // //stats_results_fout.width(10); // //stats_results_fout.right; // //stats_results_fout << hex << left ; // //stats_results_fout.setf ( ios::left, ios::adjustfield ); // //stats_results_fout.flags ( ios::right );//| ios::hex ); // //stats_results_fout << (unsigned long int)randomNo(4294967295) << endl; // //stats_results_fout.fill('0');//ios::left );//| ios::hex ); // //stats_results_fout << (unsigned long int)randomNo(4294967295) << endl; // //stats_results_fout.flags ( ios::internal | ios::hex ); // //stats_results_fout << (unsigned long int)randomNo(4294967295) << endl; // //stats_results_fout << hex << right ; // //stats_results_fout.setf ( ios::right, ios::adjustfield ); // //stats_results_fout << (unsigned long int)randomNo(4294967295) << endl; // for (total_particle_end_of_run=1 ; total_particle_end_of_run <= total_particle_per_roi; total_particle_end_of_run++) // //coef_rod_x = randomNo(1); // { // for (char k=0; k<10; k++) // { // unsigned int j = randomNo_int(); // //stats_results_fout.width(10); // //stats_results_fout.fill('0'); // if ( j < 16) // stats_results_fout << "0000000" << hex << j; // else if ( j < 256) // stats_results_fout << "000000" << hex << j; // else if ( j < 4096) // stats_results_fout << "00000" << hex << j; // else if ( j < 65536) // stats_results_fout << "0000" << hex << j; // else if ( j < 1048576) // stats_results_fout << "000" << hex << j; // else if ( j < 16777216) // stats_results_fout << "00" << hex << j; // else if ( j < 268435456) // stats_results_fout << "0" << hex << j; // else // stats_results_fout << hex << j; // } // stats_results_fout << endl; // } // //(unsigned long int)randomNo(4294967295) <<(unsigned long int)randomNo(4294967295) << \ // //(unsigned long int)randomNo(4294967295) <<(unsigned long int)randomNo(4294967295) << \ // //(unsigned long int)randomNo(4294967295) <<(unsigned long int)randomNo(4294967295) << \ // //(unsigned long int)randomNo(4294967295) <<(unsigned long int)randomNo(4294967295) << \ // //(unsigned long int)randomNo(4294967295) <<(unsigned long int)randomNo(4294967295) << \ // //endl; // write_to_log("end randomNo"); // stats_results_fout.close(); // write_to_log("start %d randMT 2^19937", total_particle_per_roi); // stats_results_fout.open(stats_results_file, ios::out | ios::app);// | std::ios::binary); // for (total_particle_end_of_run=1 ; total_particle_end_of_run <= total_particle_per_roi; total_particle_end_of_run++) // { // for (char k=0; k<10; k++) // { // unsigned int j = random_number_gen.randInt(); // //stats_results_fout.width(10); // //stats_results_fout.fill('0'); // if ( j < 16) // stats_results_fout << "0000000" << hex << j; // else if ( j < 256) // stats_results_fout << "000000" << hex << j; // else if ( j < 4096) // stats_results_fout << "00000" << hex << j; // else if ( j < 65536) // stats_results_fout << "0000" << hex << j; // else if ( j < 1048576) // stats_results_fout << "000" << hex << j; // else if ( j < 16777216) // stats_results_fout << "00" << hex << j; // else if ( j < 268435456) // stats_results_fout << "0" << hex << j; // else // stats_results_fout << hex << j; // } // stats_results_fout << endl; // } // //coef_rod_x = (float)random_number_gen.randMT(1); // //dist_n_disp_sq_fout << randomNo(1) << endl; // //stats_results_fout << std::ios::hex << \ // // random_number_gen.randInt() << random_number_gen.randInt() << \ // // random_number_gen.randInt() << random_number_gen.randInt() << \ // // random_number_gen.randInt() << random_number_gen.randInt() << \ // // random_number_gen.randInt() << random_number_gen.randInt() << \ // // random_number_gen.randInt() << random_number_gen.randInt() << \ // // endl; // write_to_log("end randMT"); // stats_results_fout.close(); // //period about 2^160 // //static // { // unsigned long x=123456789,y=362436069,z=521288629,w=88675123,v=886756453; // replace defaults with five random seed values in calling program // x = random_number_gen.randInt(); // y = random_number_gen.randInt(); // z = random_number_gen.randInt(); // w = random_number_gen.randInt(); // v = random_number_gen.randInt(); // //unsigned long xorshift(void) // write_to_log("start %d 2^160", total_particle_per_roi); // stats_results_fout.open(stats_results_file, ios::out | ios::app);// | std::ios::binary); // for (total_particle_end_of_run=1 ; total_particle_end_of_run <= total_particle_per_roi; total_particle_end_of_run++) // { // for (char k=0; k<10; k++) // { // unsigned long t; // t=(x^(x>>7)); // x=y; // y=z; // z=w; // w=v; // v=(v^(v<<6))^(t^(t<<13)); // //return // //coef_rod_x = (float)( (y+y+1)*v/4294967295.0 ); // //stats_results_fout << (y+y+1)*v << endl; // unsigned int j = (y+y+1)*v; // //stats_results_fout.width(10); // //stats_results_fout.fill('0'); // if ( j < 16) // stats_results_fout << "0000000" << hex << j; // else if ( j < 256) // stats_results_fout << "000000" << hex << j; // else if ( j < 4096) // stats_results_fout << "00000" << hex << j; // else if ( j < 65536) // stats_results_fout << "0000" << hex << j; // else if ( j < 1048576) // stats_results_fout << "000" << hex << j; // else if ( j < 16777216) // stats_results_fout << "00" << hex << j; // else if ( j < 268435456) // stats_results_fout << "0" << hex << j; // else // stats_results_fout << hex << j; // } // stats_results_fout << endl; // } // write_to_log("end 2^160"); // stats_results_fout.close(); // } // //period 2^8222 // { // static unsigned long Q[256],c=362436; /* choose random initial c<809430660 and 256 random 32-bit integers for Q[] */ // c = random_number_gen.randInt(809430660/2) + random_number_gen.randInt(809430660/2-1); //-1 for < not <= // for (file_number=0 ; file_number <= 255; file_number++) // Q[file_number] = random_number_gen.randInt(); // //unsigned long MWC256(void) // write_to_log("start %d 2^8222", total_particle_per_roi); // stats_results_fout.open(stats_results_file, ios::out | ios::app);// | std::ios::binary); // for (total_particle_end_of_run=1 ; total_particle_end_of_run <= total_particle_per_roi; total_particle_end_of_run++) // { // for (char k=0; k<10; k++) // { // unsigned long long t, a=809430660LL; // static unsigned char i=255; // t=a*Q[++i]+c; // c=(t>>32); // //return // Q[i]=t; // //coef_rod_x = (float)( t/4294967295.0 ); // //stats_results_fout << t << endl; // unsigned int j = t; // //stats_results_fout.width(10); // //stats_results_fout.fill('0'); // if ( j < 16) // stats_results_fout << "0000000" << hex << j; // else if ( j < 256) // stats_results_fout << "000000" << hex << j; // else if ( j < 4096) // stats_results_fout << "00000" << hex << j; // else if ( j < 65536) // stats_results_fout << "0000" << hex << j; // else if ( j < 1048576) // stats_results_fout << "000" << hex << j; // else if ( j < 16777216) // stats_results_fout << "00" << hex << j; // else if ( j < 268435456) // stats_results_fout << "0" << hex << j; // else // stats_results_fout << hex << j; // } // stats_results_fout << endl; // } // write_to_log("end 2^8222"); // stats_results_fout.close(); // } // //period 2^131104 // //initial seed c<809430660 and the 4096 elements (0 <= x_i < 2^32-1) in the static array Q[4096] should be assigned values before calls to the function CMWC4096(), // //otherwise the first few thousand returned values will be zeros. // //static // { // unsigned long Q[4096], c=362436; // c = random_number_gen.randInt(809430660/2) + random_number_gen.randInt(809430660/2-1); //-1 for < not <= // for (file_number=0 ; file_number < 4096; file_number++) // Q[file_number] = random_number_gen.randInt(); // //unsigned long CMWC4096(void) // write_to_log("start %d 2^131104", total_particle_per_roi); // stats_results_fout.open(stats_results_file, ios::out | ios::app); // for (total_particle_end_of_run=1 ; total_particle_end_of_run <= total_particle_per_roi; total_particle_end_of_run++) // { // for (char k=0; k<10; k++) // { // unsigned long long t, a=18782LL, b=4294967295LL; // static unsigned long i=4095; // unsigned long r=b-1;//, x; // i=(i+1)&4095; // t=a*Q[i]+c; // c=(t>>32); // t=(t&b)+c; // if(t>r) // { // c++; // t=t-b; // } // //return // Q[i]=r-t; // //coef_rod_x = (float)( Q[i]/4294967295.0 ); // //stats_results_fout << Q[i] << endl; // unsigned int j = Q[i]; // //stats_results_fout.width(10); // //stats_results_fout.fill('0'); // if ( j < 16) // stats_results_fout << "0000000" << hex << j; // else if ( j < 256) // stats_results_fout << "000000" << hex << j; // else if ( j < 4096) // stats_results_fout << "00000" << hex << j; // else if ( j < 65536) // stats_results_fout << "0000" << hex << j; // else if ( j < 1048576) // stats_results_fout << "000" << hex << j; // else if ( j < 16777216) // stats_results_fout << "00" << hex << j; // else if ( j < 268435456) // stats_results_fout << "0" << hex << j; // else // stats_results_fout << hex << j; // } // stats_results_fout << endl; // } // write_to_log("end 2^131104"); // stats_results_fout.close(); // } // //break; // //HCRYPTPROV hCryptProv; // //BYTE pbData[16]; // //bool CryptGenRandom_result; // //CryptGenRandom_result = CryptGenRandom(hCryptProv, 8, pbData) ; //} if (dlg->stop_request) goto delvar; //user aborted run //header { //stats_results_fout << "sample, " << input_base_name << endl; stats_results_fout << "Input file selected, " << input_full_name <<endl; stats_results_fout << "Output, " << output_full_name <<endl; stats_results_fout << "rod, " << rod << " , PIXEL_SIZE, " << PIXEL_SIZE << \ " , TRUE_RANDOM_WALK (no FPT)?, " << dlg->CHECK_TRUE_RANDOM_WALK.GetCheck() << \ " , FPT_boundary_layer_coef, " << FPT_boundary_layer_coef << endl; stats_results_fout << "sub_ROI_is_sphere, " << sub_ROI_is_sphere << \ " , SIZE_SUB_CUBES(radius), " << dlg->field_rw_SIZE_SUB_CUBES_float/2 << \ " , DISTANCE_FROM_CENTER, " << dlg->field_rw_DISTANCE_FROM_CENTER_MICRONS_float << ", (microns), " << endl; stats_results_fout << "mirror_coordinates, " << mirror_coordinates << \ " ,infinite, " << infinite_mirroring << " ,MAX_DISPLACEMENT_DIAG, " << (float)dlg->field_rw_MAX_DISPLACEMENT_DIAG_COEFF_float << \ " ,MAX_DISTANCE, " << (float)dlg->field_rw_MAX_DISTANCE_float << endl; stats_results_fout << "random_starting_point, " << random_starting_point << \ " ,MIN_DIST_EDGES_PIXEL, " << (float)dlg->field_rw_MIN_DIST_EDGES_PIXEL_float << endl; stats_results_fout << "check_cv, " << check_cv << \ " ,cv_num_tracers, " << cv_num_tracers << " ,cv_criteria, " << cv_criteria << " ,cv_percent_new_data, " << cv_percent_new_data << endl; stats_results_fout << "max # particles, " << total_particle_per_roi << \ " , dataCollectionPoint_increment, " << dataCollectionPoint_increment << " ,max_jump_size, " << max_jump_size << ", (microns), " << endl; stats_results_fout << "lambda, " << lambda_micron << " , lambdaForFiber, " << lambdaForFiber_micron << ", (microns), " << endl; } { coef_rod_x = (float)0; coef_rod_y = (float)0; //used to calc the average values (rod=0 case) coef_rod_z = (float)0; num_coeff = 0; coef_rod_sex = (float)0; coef_rod_sey = (float)0; coef_rod_sez = (float)0; coef_rod0_x = (float)0; coef_rod0_y = (float)0; coef_rod0_z = (float)0; } //simulation //setup the general progress bar parameters short int total_particle_for_progressbar = 1000; //scale for the progress bar displayed in dialog dlg->pbar_simulation_progress.SetRange(0, total_particle_for_progressbar); dlg->pbar_simulation_progress.SetStep(1); for (file_number=1; file_number<=total_file_number; file_number++) //do total_file_number times { //setup the progress bar parameters dlg->pbar_simulation_progress.SetPos(0); sprintf(win_txt, "Simulation: rep %d/%d (%d particles)", file_number, total_file_number, total_particle_per_roi); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); short int x_start, y_start, z_start; //voxel coord -> int {//to avoid "Warning C4701: potentially uninitialized local variable 'x_start' used" compiler warning message x_start = 0; y_start = 0; z_start = 0; } int particle; int count_tracer_start_hit_pore, count_tracer_start_hit; //used to get porosity from where (coord) the tracers land at their creation count_tracer_start_hit_pore = 0; count_tracer_start_hit = 0; num_jumps_pore_fpt = 0; num_jumps_pore_nonfpt = 0; num_jumps_fiber_fpt = 0; num_jumps_fiber_nonfpt = 0; //distance&displacement^2 values file char file_out[1024]; sprintf(file_out, "%s.simult_diff_%d.csv", pre_file_out, file_number); char* results_file; results_file = file_out; ofstream dist_n_disp_sq_fout; dist_n_disp_sq_fout.open(results_file, ios::out);// | ios::app); if(!dist_n_disp_sq_fout.is_open()) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving results: %s", file_out); write_to_log(text); return; } //file header dist_n_disp_sq_fout << "sample, " << file_out << endl; dist_n_disp_sq_fout << "Input file selected, " << input_full_name <<endl; dist_n_disp_sq_fout << "Output, " << output_full_name <<endl; dist_n_disp_sq_fout << "sim#, " << file_number << ", of, "<< total_file_number << endl; dist_n_disp_sq_fout << "rod, " << rod << ", max # particles, " << total_particle_per_roi << endl; dist_n_disp_sq_fout << "distance, x*x, y*y, z*z, x, y, z" << endl;// << endl; char stats_results_per_file_number_file_out[1024]; char* stats_results_per_file_number_file; ofstream stats_results_per_file_number_fout; if (rod == 0) { //stats results file sprintf(stats_results_per_file_number_file_out, "%s.simult_diff_results_%d.csv", pre_file_out, file_number); stats_results_per_file_number_file = stats_results_per_file_number_file_out; stats_results_per_file_number_fout.open(stats_results_per_file_number_file, ios::out | ios::app); if(!stats_results_per_file_number_fout.is_open()) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving stats results: %s", stats_results_per_file_number_fout); write_to_log(text); return; } //header stats_results_per_file_number_fout << "sample, " << file_out << endl; stats_results_per_file_number_fout << "rod, " << rod << ", max # particles, " << total_particle_per_roi << endl; stats_results_per_file_number_fout << "sim#, " << file_number << ", of, "<< total_file_number << endl; } //seed the pseudorandom-number generator before calling rand. //shouldnt be seeded again for each file,moved to top of program //srand( ( (unsigned)(time(NULL)*(double)(input_cube->wZ%211)*(double)(input_cube->wY%277)/(double)(input_cube->wX%293)) )% RAND_MAX ); //211, 277 & 293 prime //for debug, test RandomNum generator //for (file_number=1 ; file_number <= 65000; file_number++) // dist_n_disp_sq_fout << randomNo(1) << endl; //MTRand(); //random STARTING (center) POINT for ROI (within ds) unsigned short int wmini; float wmini_float=0; //find smallest cube dimension { if (input_cube->wZ <= input_cube->wX) { if (input_cube->wZ <= input_cube->wY) wmini = input_cube->wZ; else wmini = input_cube->wY; } else { if (input_cube->wX <= input_cube->wY) wmini = input_cube->wX; else wmini = input_cube->wY; } } if (!random_starting_point) { x_start = (short int)((float)wmini/2 + (float)random_number_gen.randMT((float)input_cube->wX - wmini)); // y_start = (short int)((float)wmini/2 + (float)random_number_gen.randMT((float)input_cube->wY - wmini)); // this is the center of the sub ROI by a z_start = (short int)((float)wmini/2 + (float)random_number_gen.randMT((float)input_cube->wZ - wmini)); // cube-type selection (not sphere) // if wZ is NOT min size !!PROBLEM!! (->rotate cube?) -> fixed w/ wmini // problem2? if input_cube->wZ is not even: for dimensions in x & y (shorter) - fixed by (z_start + (float)input_cube->wZ/2) below if (!mirror_coordinates) { xROImin = (float)(short int)(x_start - (float)wmini/2); xROImax = (float)(short int)(x_start + (float)wmini/2); yROImin = (float)(short int)(y_start - (float)wmini/2); //not to calc them everytime yROImax = (float)(short int)(y_start + (float)wmini/2); zROImin = (float)(short int)(z_start - (float)wmini/2); zROImax = (float)(short int)(z_start + (float)wmini/2); } else {//only for sphere type of subvolume float temp_radius; temp_radius = dlg->field_rw_SIZE_SUB_CUBES_float/(float)2; xROImin = (float)(short int)(x_start - temp_radius); xROImax = (float)(short int)(x_start + temp_radius); yROImin = (float)(short int)(y_start - temp_radius); //not to calc them everytime yROImax = (float)(short int)(y_start + temp_radius); zROImin = (float)(short int)(z_start - temp_radius); zROImax = (float)(short int)(z_start + temp_radius); } } else {//set wmini as the min(user decided min distance to edges, and smallest cube length/3) wmini_float = min((float)wmini/(float)3, (float)dlg->field_rw_MIN_DIST_EDGES_PIXEL_float); } //calculate ROI porosity float ROI_porosity=0; //unsigned char permut=0; //attempt to get rid of problems generated by bad VS2005 internal pseudo random gen //if (!random_starting_point) { unsigned long long ROI_porosity_fiber_pores_long = 0; unsigned long long ROI_porosity_total_pixels_long = 0; if (!mirror_coordinates & !random_starting_point) { if (!sub_ROI_is_sphere) {//sub-ROI = cube case unsigned short int x, y, z; for (z = (short int)zROImin; z < (short int)zROImax; z++) for (y = (short int)yROImin; y < (short int)yROImax; y++) for (x = (short int)xROImin; x < (short int)xROImax; x++) { ROI_porosity_total_pixels_long ++; if ( !input_cube->get_spot(z, y, x) ) //no need to modify the xyz, they cannot be mirrored in the sub-roi = cube case ROI_porosity_fiber_pores_long ++; } } else {// sphere type of sub-volume int x, y, z; for (z = (short int)zROImin; z < (short int)zROImax; z++) for (y = (short int)yROImin; y < (short int)yROImax; y++) for (x = (short int)xROImin; x < (short int)xROImax; x++) if ( pow(x-(float)x_start, 2) \ + pow(y-(float)y_start, 2) \ + pow(z-(float)z_start, 2) \ < sq_MAX_ROI_RADIUS ) { ROI_porosity_total_pixels_long ++; if ( !input_cube->get_spot( \ (coord_trans_int((z), input_cube->wZ)), \ (coord_trans_int((y), input_cube->wY)), \ (coord_trans_int((x), input_cube->wX)) ) ) ROI_porosity_fiber_pores_long ++; } } } else //MIRROR CASE or random starting point -> consider cube case/ disregard sphere shape {// sphere type of sub-volume unsigned short int x, y, z; for (z = 0; z < input_cube->wZ; z++) for (y = 0; y < input_cube->wY; y++) for (x = 0; x < input_cube->wX; x++) { ROI_porosity_total_pixels_long ++; if ( !input_cube->get_spot(z, y, x) ) ROI_porosity_fiber_pores_long ++; } } ROI_porosity = (float)ROI_porosity_fiber_pores_long / (float)ROI_porosity_total_pixels_long; if ( (ROI_porosity == 0) && (rod == 0) ) { write_to_log("porosity=0 AND rod=0 -> skipping to next random ROI/sub-cube"); stats_results_per_file_number_fout << "porosity=0 AND rod=0 -> skipping to next random ROI/sub-cube" << endl; if (rod == 0) stats_results_per_file_number_fout.close(); //close stats values file continue; } //ROI_porosity /= pow((float)wmini, 3); total_ROI_porosity += ROI_porosity; //write_to_log("true_ROI_porosity= %f", ROI_porosity); stats_results_fout << "***sim#, " << file_number << ", true_ROI_porosity, " << ROI_porosity << endl; stats_results_fout << "subROIcenter, " << x_start << ", " << y_start << ", " << z_start << endl; if ( (rod == 0) )//&& dlg->CHECK_SAVE_DETAILED_TRACE.GetCheck() ) { stats_results_per_file_number_fout << "***sim#, " << file_number << ", ROI_porosity, " << ROI_porosity << endl; stats_results_per_file_number_fout << "xyz_start, " << x_start << ", " << y_start << ", " << z_start << endl; } } //save tracer path - for debug DS_1b_cube * ds_cube_tracer_path; //for (particle=1; particle<=total_particle_per_roi; particle++) //do total_particle_per_roi times particle=1; do { //if user requested termination of simulation //if (dlg->stop_request) // goto delvar; sprintf(win_txt, "Simulation: rep %d/%d (%d)", file_number, total_file_number, particle); dlg->txt_field_progress.SetWindowText(win_txt); float x, y, z; //tracer coords -> float float xo, yo, zo; //initial tracer coord float x2, y2, z2; x2=0; y2=0; z2=0; //float xf, yf, zf; //final coord after exiting ROI float r; //radius/euclidian distance till interface //short int xo, yo, zo; //initial tracer coord -> int float dataCollectionPoint; dataCollectionPoint = (float)dataCollectionPoint_increment; distance = 0; //searchStartingPoint (); if (!random_starting_point) //random STARTING POINT for tracer, but close to center of sub volume (not real random case with infinite mirroring) { if (!sub_ROI_is_sphere) {//this is a cube type of sub volume r = dlg->field_rw_DISTANCE_FROM_CENTER_float; x = xo = ( x_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); y = yo = ( y_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); // randomNo(0.0001) to avoid int numbers z = zo = ( z_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); //which can cause some issues further down } else {//sphere type of sub volume float alpha, beta, gamma, mod; do //hypercube rejection method { alpha = 1-2*(float)random_number_gen.randExc(); beta = 1-2*(float)random_number_gen.randExc(); gamma = 1-2*(float)random_number_gen.randExc(); mod = (alpha*alpha + beta*beta + gamma*gamma); } while (mod < 0.0001 || mod > 1); mod = sqrt(mod); r = (float)random_number_gen.randMT(dlg->field_rw_DISTANCE_FROM_CENTER_float); //max radius x = xo = ( x_start - r*(alpha/mod) ) + (float)random_number_gen.randMT((float)0.001); y = yo = ( y_start - r* (beta/mod) ) + (float)random_number_gen.randMT((float)0.001); z = zo = ( z_start - r*(gamma/mod) ) + (float)random_number_gen.randMT((float)0.001); } {//update counters for ROI porosity based on where tracers start count_tracer_start_hit_pore += (!input_cube->get_spot(\ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ); count_tracer_start_hit ++; } if (rod == 0) {//issue: what if there is no pore voxel within selected ROI and rod==0? -> detected when calculating ROI_porosity & avoided // but ok in 'simultaneous' diffusion to have the tracer start in fiber space unsigned short int counter = 0; //to avoid infinite loop while ( input_cube->get_spot(\ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ) { if (!sub_ROI_is_sphere) {//this is a cube type of sub volume x = xo = ( x_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); y = yo = ( y_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); z = zo = ( z_start - r*(1-2*(float)random_number_gen.randMT(1)) ) + (float)random_number_gen.randMT((float)0.001); }// loop to ensure starting point is in pore space else {//hypercube rejection method -> sphere type of sub volume counter++; float alpha, beta, gamma, mod; do { alpha = 1-2*(float)random_number_gen.randExc(); beta = 1-2*(float)random_number_gen.randExc(); gamma = 1-2*(float)random_number_gen.randExc(); mod = (alpha*alpha + beta*beta + gamma*gamma); } while (mod < 0.0001 || mod > 1); mod = sqrt(mod); //problem: infinite loop if rod=0 and dlg->field_rw_DISTANCE_FROM_CENTER_float=0 // and the center pixel is fiber -> need to increase the radius size if (counter > 1000) { //r = randomNo( .5 * (float)dlg->field_rw_SIZE_SUB_CUBES_float ); //max possible value counter = 1000; //to avoid overflow counter variable } //else { if (dlg->field_rw_DISTANCE_FROM_CENTER_float == 0) //r = randomNo( .1 * dlg->field_rw_SIZE_SUB_CUBES_float ); //fixed @ 10% of subROI size -> might not be enough!! r = (float)random_number_gen.randMT( (float)counter/(float)2000 * dlg->field_rw_SIZE_SUB_CUBES_float ); //lin var with counter //r = randomNo( 2*pow((float)counter/(float)2000, 2) * dlg->field_rw_SIZE_SUB_CUBES_float ); //non-lin var with counter else r = (float)random_number_gen.randMT(dlg->field_rw_DISTANCE_FROM_CENTER_float); //max radius } x = xo = ( x_start - r*(alpha/mod) ) + (float)random_number_gen.randMT((float)0.001); y = yo = ( y_start - r* (beta/mod) ) + (float)random_number_gen.randMT((float)0.001); z = zo = ( z_start - r*(gamma/mod) ) + (float)random_number_gen.randMT((float)0.001); } {//update counters for ROI porosity based on where tracers start count_tracer_start_hit ++; //need to take into account every time new coord are generated //count_tracer_start_hit_pore++ //no need, has to be a pore to exit do loop } } } } else //true random STARTING POINT, used with infinite mirroring {//wmini_float is already the min(min dist to edges & wmini/3) x = xo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wX - 2*wmini_float)); // y = yo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wY - 2*wmini_float)); // this is selected in the sub ROI by a z = zo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wZ - 2*wmini_float)); // cube-type selection (not sphere) {//update counters for ROI porosity based on where tracers start count_tracer_start_hit_pore += (!input_cube->get_spot(\ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ); count_tracer_start_hit ++; } if (rod == 0) {//issue: what if there is no pore voxel within selected ROI and rod==0? //-> could happen if all pores are close to edges (unlikely) - because of: wmini_float = min(... above //this would create an infinite loop here. // but ok in 'simultaneous' diffusion to have the tracer start in fiber space while ( input_cube->get_spot(\ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ) {//while in fiber x = xo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wX - 2*wmini_float)); // y = yo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wY - 2*wmini_float)); // this is selected in the sub ROI by a z = zo = (wmini_float + (float)random_number_gen.randMT((float)input_cube->wZ - 2*wmini_float)); // cube-type selection (not sphere) //update counters for ROI porosity based on where tracers start count_tracer_start_hit ++; //need to take into account every time new coord are generated } } x_start = (short int)(xo + .5); y_start = (short int)(yo + .5); z_start = (short int)(zo + .5); } if ( dlg->CHECK_SAVE_DETAILED_TRACE.GetCheck() ) { dist_n_disp_sq_fout << "**particle#, " << particle; dist_n_disp_sq_fout << ", @, coord, " << xo << ", " << yo << ", " << zo << endl; //save tracer path - for debug //DS_1b_cube * ds_cube_tracer_path; ds_cube_tracer_path = new DS_1b_cube( input_cube->wZ, input_cube->wY, input_cube->wX ); // Initial state of tracer path: empty ds_cube_tracer_path->clear(); //save tracer position into path map ds_cube_tracer_path->set_spot_1( (unsigned short int)z, (unsigned short int)y, (unsigned short int)x ); } while (distance < dlg->field_rw_MAX_DISTANCE_float) { if ( !input_cube->get_spot(\ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX))) ) {//tracer in pore space bool fpt_flag = false; if ( !dlg->CHECK_TRUE_RANDOM_WALK.GetCheck() ) //FPT allowed { //findNearestFiber (); // gets r (max) at that point //temp buffer for min radius determination //float * radius_values = new float [27]; //26 3D-neighbours + center voxel //kind of convolution char i,j,k; //unsigned char ind_val = 0; //float minimum; r = 1 + (float)distMap->datac( \ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ; //min radius cannot be bigger than @ location in map + 1 for (i=-1; i<2; i++) for (j=-1; j<2; j++) for (k=-1; k<2; k++) { if ( !mirror_coordinates && ( ((int)x+(int)i < 0) || ((int)x+(int)i > input_cube->wX-1) || ((int)y+(int)j < 0) || ((int)y+(int)j > input_cube->wY-1) || ((int)z+(int)k < 0) || ((int)z+(int)k > input_cube->wZ-1)) ) {//coord out of input_cube -> do nothing //radius_values[ind_val] = (float)distMap->datac((unsigned short int)z, (unsigned short int)y, (unsigned short int)x) + 1; //ind_val++; //min radius cannot be bigger than @ location in map + 1 } else {//valid coord float distance_temp; distance_temp = (float)distMap->datac(\ coord_trans_int((int)floor(z+k), input_cube->wZ), \ coord_trans_int((int)floor(y+j), input_cube->wY), \ coord_trans_int((int)floor(x+i), input_cube->wX) ) \ * !input_cube->get_spot(\ coord_trans_int((int)floor(z+k), input_cube->wZ), \ coord_trans_int((int)floor(y+j), input_cube->wY), \ coord_trans_int((int)floor(x+i), input_cube->wX) ); //radius = 0 if voxel is fiber float x_fracpart=0, y_fracpart=0, z_fracpart=0; //fractional parts of coords for DISTANCE to convolution voxel //float intpart; //x_fracpart = abs(modf (x , &intpart)); //fractional part of x, with the same sign //if (i>0) // x_fracpart = 1 - x_fracpart; //y_fracpart = modf (y , &intpart); //if (j>0) // y_fracpart = 1 - y_fracpart; //z_fracpart = modf (z , &intpart); //if (k>0) // z_fracpart = 1 - z_fracpart; // = fracpart if x>0, =1-fracpart if x<0 x_fracpart = x - (float)(int)floor(x); if (i>0) x_fracpart = (float)1 - x_fracpart; y_fracpart = y - (float)(int)floor(y); if (j>0) y_fracpart = (float)1 - y_fracpart; z_fracpart = z - (float)(int)floor(z); if (k>0) z_fracpart = (float)1 - z_fracpart; //if (i>0) //this is wrong // x_fracpart = (float)1 - x_fracpart; //else // x_fracpart = x - (float)(int)floor(x); //if (j>0) // y_fracpart = (float)1 - y_fracpart; //else // y_fracpart = y - (float)(int)floor(y); //if (k>0) // z_fracpart = (float)1 - z_fracpart; //else // z_fracpart = z - (float)(int)floor(z); //if (i>0) // x_fracpart = 1-(coord_trans_float(x, input_cube->wX)-(int)coord_trans_float(x, input_cube->wX)); //else // x_fracpart = (coord_trans_float(x, input_cube->wX)-(int)coord_trans_float(x, input_cube->wX)); //if (j>0) // y_fracpart = 1-(coord_trans_float(y, input_cube->wY)-(int)coord_trans_float(y, input_cube->wY)); //else // y_fracpart = (coord_trans_float(y, input_cube->wY)-(int)coord_trans_float(y, input_cube->wY)); //if (k>0) // z_fracpart = 1-(coord_trans_float(z, input_cube->wZ)-(int)coord_trans_float(z, input_cube->wZ)); //else // z_fracpart = coord_trans_float(z, input_cube->wZ)-(int)coord_trans_float(z, input_cube->wZ); if (distance_temp == 0) //simple Pythagore, only in the direction of the target pixel { //radius_values[ind_val] = sqrt( pow(x_fracpart,2)+pow(y_fracpart,2)+pow(z_fracpart,2) ); if ( (i==0) + (j==0) + (k==0) == 3) continue; else distance_temp = sqrt( (i!=0)*pow(x_fracpart,2) \ + (j!=0)*pow(y_fracpart,2) \ + (k!=0)*pow(z_fracpart,2) ); //if (r > distance_temp) // r = distance_temp; } else switch ((i==0) + (j==0) + (k==0)) { case 2: if (i) //radius_values[ind_val] = x_fracpart; distance_temp += x_fracpart; else if (j) //radius_values[ind_val] = y_fracpart; distance_temp += y_fracpart; else if (k) //radius_values[ind_val] = z_fracpart; distance_temp += z_fracpart; //radius_values[ind_val] += distance_temp; //if (r > distance_temp) // r = distance_temp; break; case 1: if (i==0) { //determine coord closest to convolution voxel, then calc min possible radius (Pythagore) if ( y_fracpart < z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(y_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(y_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); else //z closest //radius_values[ind_val] = sqrt( pow(z_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(z_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); } else if (j==0) { if ( x_fracpart < z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(x_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); else //z closest //radius_values[ind_val] = sqrt( pow(z_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(z_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); } else if (k==0) { if ( x_fracpart < y_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(x_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); else //y closest //radius_values[ind_val] = sqrt( pow(y_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(y_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); } //if (r > distance_temp) // r = distance_temp; break; case 0: //determine coord closest to convolution voxel, then calc if (distance_temp == 1) { if ( x_fracpart <= y_fracpart ) //x closest { if ( x_fracpart <= z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart+1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); } else //y closest { if ( y_fracpart <= z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+1, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+1, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); } } else { if ( x_fracpart <= y_fracpart ) //x closest { if ( x_fracpart <= z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+distance_temp-1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart+distance_temp-1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); } else //y closest { if ( y_fracpart <= z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+distance_temp-1, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+distance_temp-1, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); } } //if (r > distance_temp) // r = distance_temp; break; default: //case 3: //center voxel, ind_val=13? -> discard value //radius_values[ind_val] = (float)distMap->datac((unsigned short int)z, (unsigned short int)y, (unsigned short int)x) + 1; //min radius cannot be bigger than @ location in map + 1 //radius_values[ind_val] = 9999999; //ignore center voxel by setting large value -> not min distance_temp = r; break; } //ind_val++; if (r > distance_temp) r = distance_temp; //get minimum value for r } } //r = mini(radius_values, 27); //get min of the min radius //delete [] radius_values; r = min(r, max_jump_size_lambda_pix); //added for debug: for empty samples, need to limit size of jump in FPT? if (r <= FPT_boundary_layer_lambda_pix) // if radius < 5 * mean free pass // 5*lambda~.7um //fpt_flag = true; //while (r > FPT_boundary_layer_lambda_pix) { //do //r = -lambda_pix*log( 1-(float)random_number_gen.randMT(1) ); // 1.0000001 to avoid #IND after log <- no need r = -lambda_pix*(float)log( 1-random_number_gen.randExc() ); // randExc <- needed! //fpt_flag = false; //dist_n_disp_sq_fout << r << endl; //while ( r>=1 ); // -> to avoid problems when finding interface - limit size of random jump // but problems can still happen (need to check all voxels traversed by tracer path between old and new coords) } else { fpt_flag = true; //r -= FPT_boundary_layer_lambda_pix_4_fifth; //to land at min 4/5*lambda away from the interface //commented out for small likelyhood } //if (r!=r) //for debug, #IND tracking quiet_NaN // dist_n_disp_sq_fout << r << "," ; //debug } else //forced TRUE_RANDOM_WALK mode - no FPT allowed { r = -lambda_pix*(float)log( 1-random_number_gen.randExc() ); //r = min(r, max_jump_size*lambda_pix); // to limit size of jump } //jumpinPore (); // random jump (now that min r is known) {//tracer in pore float alpha, beta, gamma, mod; //, sign; float dx, dy, dz; do //hypercube rejection method { alpha = 1-2*(float)random_number_gen.randExc();//randMT(1); beta = 1-2*(float)random_number_gen.randExc();//randMT(1); gamma = 1-2*(float)random_number_gen.randExc();//randMT(1); mod = (alpha*alpha + beta*beta + gamma*gamma); } while (mod < 0.0001 || mod > 1); mod = sqrt(mod); //switch (permut%3) //{//attempt to get rid of non-randomness //case 0: dx = r*(alpha/mod); dy = r*(beta/mod); //jump coord dz = r*(gamma/mod); // break; //case 1: // dz = r*(alpha/mod); // dx = r*(beta/mod); // dy = r*(gamma/mod); // break; //default: //case 2: // dy = r*(alpha/mod); // dz = r*(beta/mod); // dx = r*(gamma/mod); // break; //} //permut ++; x2 = dx + x; y2 = dy + y; //hence new coordinates z2 = dz + z; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } //if ( (!fpt_flag) && (input_cube->get_spot( if ( (!fpt_flag) && (input_cube->get_spot( \ coord_trans_int((int)floor(z2), input_cube->wZ), \ coord_trans_int((int)floor(y2), input_cube->wY), \ coord_trans_int((int)floor(x2), input_cube->wX) ) ) ) //tracer (new coord) belongs to fiber space {//should only happen in *real/pure* random walk (fpt_flag=0) case and when (rod != 0) //need to check if we skipped some FIBER pixels in between, and locate the interface //check for every voxel on the path between origin (x)-pore and end (x2)-fiber voxels //- stop as soon as we meet a FIBER voxel, to get the corresponding interface coord float r_index_final; float xi, yi, zi; //interface coord (should be int for at least one of them) { r_index_final = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x2 > x ) {//********* for (x_index = (int)ceil(x); x_index <= (int)floor(x2); x_index++) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); // y = old + ratio * (new - old) z_index = (int)floor(z + r_index*(z2-z)); //check if coord valid (belongs between points X and X2 - rounded to closest int) if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if fiber r_index_final = r_index; break; } } } else { for (x_index = (int)floor(x); x_index >= (int)ceil(x2); x_index--) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); z_index = (int)floor(z + r_index*(z2-z)); if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final = r_index; break; } } } //for Y if ( y2 > y ) { for (y_index = (int)ceil(y); y_index <= (int)floor(y2); y_index++) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (y_index = (int)floor(y); y_index >= (int)ceil(y2); y_index--) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(y_index+1) - y) / (y2 - y)); break; } } } //for Z if ( z2 > z ) { for (z_index = (int)ceil(z); z_index <= (int)floor(z2); z_index++) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (z_index = (int)floor(z); z_index >= (int)ceil(z2); z_index--) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(z_index+1) - z) / (z2 - z)); break; } } } //if (r_index_final != (float)1) //re calc the coord to be at the interface - coord(i) xi = (x + r_index_final*(x2-x)); yi = (y + r_index_final*(y2-y)); zi = (z + r_index_final*(z2-z)); } //now decide if reflection or enter fiber {//original tracer coords (in pore space) before random jump float x3, y3, z3; //new coord after taking into account the interface coords float x4=0, y4=0, z4=0; //new coord after taking into account the possible jump over pore if ((float)random_number_gen.randExc() < sqrt_rod) //should not be <= otherwise could enter fiber in rod=0 case (with probability = 0!) but taken care of by def of sqrt_rod = -1 when rod =0 { // pore -> enter fiber space x3 = xi + sqrt_rod*(x2-xi); y3 = yi + sqrt_rod*(y2-yi); //going from high mobility media to low mobility media z3 = zi + sqrt_rod*(z2-zi); //if (rod>1) { //need to check we didnt jump over some pore voxels (mostly rod>1 case - but could happen anyway) //check for every voxel on the path between origin (xi)-fiber and end (x3)-fiber or pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final2; float xi2, yi2, zi2; //interface coord (should be int for at least one of them) { r_index_final2 = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x3 > xi ) { for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) //start @ xi+1, finish @ x3 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) z_index = (int)floor(zi + r_index*(z3-zi)); //check if coord valid in input_cube if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if pore r_index_final2 = r_index; break; } } } } else { for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) //start @ xi, finish @ x3+1 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final2 = r_index; break; } } } } //for Y if ( y3 > yi ) { for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } //for Z if ( z3 > zi ) { for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } ////for X //if ( x3 > xi ) //{ // for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) //start @ xi+1, finish @ x3 // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) // z_index = (int)floor(zi + r_index*(z3-zi)); // //check if coord valid in input_cube // if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && // (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // {//if pore // r_index_final2 = r_index; // break; // } // } //} //else //{ // for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) //start @ xi, finish @ x3+1 // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && // (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index-1, input_cube->wX) ) ) // { // r_index_final2 = r_index; // break; // } // } //} ////for Y //if ( y3 > yi ) //{ // for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && // (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final2 = min(r_index_final2, r_index); // break; // } // } //} //else //{ // for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && // (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index-1, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final2 = min(r_index_final2, r_index); // break; // } // } //} ////for Z //if ( z3 > zi ) //{ // for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && // (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final2 = min(r_index_final2, r_index); // break; // } // } //} //else //{ // for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && // (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) // if (!input_cube->get_spot( \ // coord_trans_int(z_index-1, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final2 = min(r_index_final2, r_index); // break; // } // } //} if (r_index_final2 != (float)1) { //re calc the coord to be at the 2nd interface - coord(i2) xi2 = (xi + r_index_final2*(x3-xi)); yi2 = (yi + r_index_final2*(y3-yi)); zi2 = (zi + r_index_final2*(z3-zi)); //then jump close to interface #2 minus a small distance (lambda?) r_index = sqrt( pow(xi2-xi,2) + pow(yi2-yi,2) + pow(zi2-zi,2) ); //distance between first interface and second one if (r_index > two_FPT_boundary_layer_lambda_pix) //after jump should be at least lambda away from r_index -= FPT_boundary_layer_lambda_pix; else r_index *= (float)0.55; //to be slightly closer to the 2nd interface, yet as far as possible to respect the boundary_layer_lambda distance //r_index = min(r_index, ); //r_index = min(r_index, ); x4 = xi2 + r_index * (xi - xi2); y4 = yi2 + r_index * (yi - yi2); z4 = zi2 + r_index * (zi - zi2); } else {//no jump over pore, coord are valid, keep them x4 = x3; y4 = y3; z4 = z3; } } } //else //rod<1 //{//no jump over pore, coord are valid, keep them // x4 = x3; // y4 = y3; // z4 = z3; //} //if ( !input_cube->get_spot( \ // (unsigned short int)coord_trans_int((int)floor(z3), input_cube->wZ), \ // (unsigned short int)coord_trans_int((int)floor(y3), input_cube->wY), \ // (unsigned short int)coord_trans_int((int)floor(x3), input_cube->wX) ) ) //{//re-pool the distribution // #if RW_DO_DEBUG // write_to_log("error: should enter fiber space, ->pore instead, will now 'continue'"); // #endif // continue; //} //re-calculate r because different (shorter if rod<=1) than previously calculated (new=x4, previous before jump =x) r = sqrt( pow(x4-x, 2)+pow(y4-y, 2)+pow(z4-z, 2)); //calc distance //distance += r; //in pixels #if RW_DO_DEBUG if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug {//should not happen, not in FPT, unless wrong calc of r (too large) write_to_log("error FPT -1-"); } #endif distance += (float)PIXEL_SIZE*r; //in microns //should separate the part within pore space and the part within fiber space? -> debatable //no need to check if coord 3 is valid, because coord 2 is -> wrong, but taken care above //save new coords as being current ones x = x4; y = y4; z = z4; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } else { // reflected back (pore -> pore) //reflectFlag = 1; //{//bounce of how much would have entered freely // x3 = xi - (x2-xi); // y3 = yi - (y2-yi); //can cause problems (large bounce) // z3 = zi - (z2-zi); //} //{//bounce half way till the interface // x3 = x + (xi-x)/2; // y3 = y + (yi-y)/2; // z3 = z + (zi-z)/2; //} {//bounce back a random distance from the interface, depending on fractional_part(r) float fractpart, intpart; fractpart = modf (r , &intpart); //if = 0 then problem! can enter fiber (with x being integer) - very unlikely - taken care of with the 'continue" below x3 = x + (xi-x)*fractpart; y3 = y + (yi-y)*fractpart; z3 = z + (zi-z)*fractpart; } if ( input_cube->get_spot( \ (unsigned short int)coord_trans_int((int)floor(z3), input_cube->wZ), \ (unsigned short int)coord_trans_int((int)floor(y3), input_cube->wY), \ (unsigned short int)coord_trans_int((int)floor(x3), input_cube->wX) ) ) {//re-pool the exp distribution #if RW_DO_DEBUG write_to_log("error reflection1 ->fiber instead of pore, will now 'continue'"); //debug #endif continue; } //re-calculate r because shorter than previously calculated (new=x3, previous before jump =x) //NOT taking into account the bouncing on the interface: //r = sqrt( pow(x3-x, 2)+pow(y3-y, 2)+pow(z3-z, 2) ); //taking into account the bouncing on the interface: //-> 150% distance between x and interface (sqrt(1.5)=1.224745) //r = (float)1.224745*sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ); //150%, bounce 1/2 way //taking into account the bouncing on the interface: r = sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ) \ + sqrt( pow(xi-x3, 2)+pow(yi-y3, 2)+pow(zi-z3, 2) ); // ||x,xi||+||xi,x3|| //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug //{//should not happen, not in FPT, unless wrong calc of r // #if RW_DO_DEBUG // write_to_log("error FPT -2-"); // #endif // //distance += pow(r, 2)/(2*lambda_pix); //in pixels // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns // distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //} //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns //new coord (x3) is valid because x2 is when bouncing back //-> no, because of the reflection formula (the first one) //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // { // //write_to_log("simultaneous diffusion: pore->fiber -- ERROR2 -- "); // //if tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} //save new coords as being current ones x = x3; y = y3; z = z3; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } } } else //tracer (new coord) belongs to pore space { //if (input_cube->get_spot( coord_trans_int((int)floor(z2), input_cube->wZ), \ // coord_trans_int((int)floor(y2), input_cube->wY), \ // coord_trans_int((int)floor(x2), input_cube->wX) )) // continue; //if fiber //for debug if ( fpt_flag ) {//keep new coordinates, tracer did NOT jump over some FIBER voxels since we are in FPT case x = x2; y = y2; z = z2; //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) //distance += pow(r, 2)/(2*lambda_pix); //in pixels //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //else // //distance += r; //in pixels // distance += PIXEL_SIZE*r; //in microns } else {//need to check the tracer didnt jump over some fiber voxels, can happen if (random(r) > (r calc for fpt) but this is a non FPT case) //need to check if we skipped some FIBER pixels in between, and locate the interface //check for every voxel on the path between origin (x)-pore and end (x2)-fiber voxels //- stop as soon as we meet a FIBER voxel, to get the corresponding interface coord float r_index_final; float xi, yi, zi; //interface coord (should be int for at least one of them) r_index_final = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x2 > x ) { for (x_index = (int)ceil(x); x_index <= (int)floor(x2); x_index++) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); // y = old + ratio * (new - old) z_index = (int)floor(z + r_index*(z2-z)); //check if coord valid ( min <= x_index < max ) if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if fiber r_index_final = r_index; //no need for min since we "break" asa it happens break; } } } else { for (x_index = (int)floor(x); x_index >= (int)ceil(x2); x_index--) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); z_index = (int)floor(z + r_index*(z2-z)); if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final = r_index;//((float)(x_index+1) - x) / (x2 - x); break; } } } //for Y if ( y2 > y ) { for (y_index = (int)ceil(y); y_index <= (int)floor(y2); y_index++) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (y_index = (int)floor(y); y_index >= (int)ceil(y2); y_index--) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(y_index+1) - y) / (y2 - y)); break; } } } //for Z if ( z2 > z ) { for (z_index = (int)ceil(z); z_index <= (int)floor(z2); z_index++) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (z_index = (int)floor(z); z_index >= (int)ceil(z2); z_index--) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(z_index+1) - z) / (z2 - z)); break; } } } if (r_index_final != (float)1) {//tracer did jump over fiber voxels, find interface #if RW_DO_DEBUG if ( fpt_flag ) //for debug write_to_log("error fpt jumped over fiber"); #endif //re calc the coord to be at the interface - coord(i) //if ( fpt_flag ) // continue; //debug xi = (x + r_index_final*(x2-x)); yi = (y + r_index_final*(y2-y)); zi = (z + r_index_final*(z2-z)); //now decide if reflection or enter fiber {//original tracer coords (in pore space) before random jump float x3, y3, z3; //new coord after taking into account the interface coords float x4=0, y4=0, z4=0; //new coord after taking into account the possible jump over pore if ((float)random_number_gen.randExc() < sqrt_rod) // pore -> enter fiber space { x3 = xi + sqrt_rod*(x2-xi); y3 = yi + sqrt_rod*(y2-yi); //going from high mobility media to low mobility media z3 = zi + sqrt_rod*(z2-zi); //if (rod>1) { //need to check we didnt jump over some pore voxels (rod>1 case - no can also happen w/ rod<=1) //check for every voxel on the path between origin (xi)-fiber and end (x3)-fiber or pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final2; float xi2, yi2, zi2; //interface coord (should be int for at least one of them) { r_index_final2 = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x3 > xi ) { for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) //start @ xi+1, finish @ x3 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) z_index = (int)floor(zi + r_index*(z3-zi)); //check if coord valid in input_cube if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if pore r_index_final2 = r_index; break; } } } } else { for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) //start @ xi, finish @ x3+1 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final2 = r_index; break; } } } } //for Y if ( y3 > yi ) { for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } //for Z if ( z3 > zi ) { for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } if (r_index_final2 != (float)1) { //re calc the coord to be at the interface - coord(i) xi2 = (xi + r_index_final2*(x3-xi)); yi2 = (yi + r_index_final2*(y3-yi)); zi2 = (zi + r_index_final2*(z3-zi)); //then jump close to interface #2 minus a small distance (lambda?) r_index = sqrt( pow(xi2-xi,2) + pow(yi2-yi,2) + pow(zi2-zi,2) ); //distance between first interface and second one if (r_index > two_FPT_boundary_layer_lambda_pix) //after jump should be at least 4/5*lambda away from interface r_index -= FPT_boundary_layer_lambda_pix; //FPT_boundary_layer_lambda_pix; else r_index *= (float)0.55; //r_index = min(r_index, ); //r_index = min(r_index, ); x4 = xi2 + r_index * (xi - xi2); y4 = yi2 + r_index * (yi - yi2); z4 = zi2 + r_index * (zi - zi2); } else {//no jump over pore, coord are valid, keep them x4 = x3; y4 = y3; z4 = z3; } } //if ( !input_cube->get_spot( \ // (unsigned short int)coord_trans_int((int)floor(z3), input_cube->wZ), \ // (unsigned short int)coord_trans_int((int)floor(y3), input_cube->wY), \ // (unsigned short int)coord_trans_int((int)floor(x3), input_cube->wX) ) ) //{//re-pool the distribution // #if RW_DO_DEBUG // write_to_log("error: should enter fiber space, ->pore instead, will now 'continue'"); // #endif // continue; //} } //else //{//no jump over pore, coord are valid, keep them // x4 = x3; // y4 = y3; // z4 = z3; //} //re-calculate r because different (shorter if rod<=1) than previously calculated (new=x4, previous before jump =x) r = sqrt( pow(x4-x, 2)+pow(y4-y, 2)+pow(z4-z, 2)); //re-calculate r because shorter (only if rod<=1) than previously calculated (new=x3, previous before jump =x) //r = sqrt( pow(x3-x, 2)+pow(y3-y, 2)+pow(z3-z, 2)); //calc distance //distance += r; //in pixels #if RW_DO_DEBUG if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug {//should not happen, not in FPT, unless wrong calc of r (too large) write_to_log("error FPT -1-"); } #endif distance += (float)PIXEL_SIZE*r; //in microns //should separate the part within pore space and the part within fiber space? //no need to check if coord 3 is valid, because coord 2 is } else // reflected back (pore -> pore) { //reflectFlag = 1; //{//bounce of how much would have entered freely // x3 = xi - (x2-xi); // y3 = yi - (y2-yi); //can cause problems (large bounce) // z3 = zi - (z2-zi); //} //{//bounce half way till the interface // x3 = x + (xi-x)/2; // y3 = y + (yi-y)/2; // z3 = z + (zi-z)/2; //} {//bounce random way till the interface, depending on fractional_part(r) float fractpart, intpart; fractpart = modf (r , &intpart); x3 = x + (xi-x)*fractpart; y3 = y + (yi-y)*fractpart; z3 = z + (zi-z)*fractpart; } if ( input_cube->get_spot( \ coord_trans_int((int)floor(z3), input_cube->wZ), \ coord_trans_int((int)floor(y3), input_cube->wY), \ coord_trans_int((int)floor(x3), input_cube->wX) ) ) {//re-pool the exp distribution #if RW_DO_DEBUG write_to_log("error reflection2 ->fiber instead of ->pore, will now 'continue'"); //debug #endif continue; } //re-calculate r because shorter than previously calculated (new=x3, previous before jump =x) //NOT taking into account the bouncing on the interface: //r = sqrt( pow(x3-x, 2)+pow(y3-y, 2)+pow(z3-z, 2) ); //taking into account the bouncing on the interface: //-> 150% distance between x and interface (sqrt(1.5)=1.224745) //r = (float)1.224745*sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ); r = sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ) \ + sqrt( pow(xi-x3, 2)+pow(yi-y3, 2)+pow(zi-z3, 2) ); // ||x,xi||+||xi,x3|| //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug //{//should not happen, not in FPT, unless wrong calc of r // #if RW_DO_DEBUG // write_to_log("error FPT -2-"); // #endif // //distance += pow(r, 2)/(2*lambda_pix); //in pixels // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns // distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //} //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns //new coord (x3) *should* be valid because x2 is (case: bounce half way till the interface) //(no, because of the reflection formula case: bounce of how much would have entered freely) //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // { // //write_to_log("simultaneous diffusion: pore->fiber -- ERROR2 -- "); // //if tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} } //else, save new coords as being current ones x = x3; y = y3; z = z3; #if RW_DO_DEBUG if ( input_cube->get_spot( (unsigned short int)z, (unsigned short int)y, (unsigned short int)x ) ) write_to_log("error fpt #7"); //debug #endif } } else {//keep new coordinates, tracer did NOT jump over some FIBER voxels x = x2; y = y2; z = z2; //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) // //distance += pow(r, 2)/(2*lambda_pix); //in pixels // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns // distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns } } } //calculateDistance (); //done above if (fpt_flag) ++num_jumps_pore_fpt; else ++num_jumps_pore_nonfpt; } } else {// tracer in fiber space -> can only be when (rod != 0) //phase transition should only happen in *real/pure* random walk case #if RW_DO_DEBUG if ( rod == (float)0 ) {// - to debug write_to_log("error rod==0 & entered fiber!, will now 'break'"); break; } #endif bool fpt_flag = false; if (!dlg->CHECK_TRUE_RANDOM_WALK.GetCheck()) //FPT allowed { //findNearestPore ( ); //temp buffer for min radius determination //float * radius_values = new float [27]; //26 3D-neighbours + center voxel //kind of convolution char i,j,k; //unsigned char ind_val = 0; r = 1 + (float)distMap->datac( \ (coord_trans_int((int)floor(z), input_cube->wZ)), \ (coord_trans_int((int)floor(y), input_cube->wY)), \ (coord_trans_int((int)floor(x), input_cube->wX)) ) ; //min radius cannot be bigger than @ location in map + 1 for (i=-1; i<2; i++) for (j=-1; j<2; j++) for (k=-1; k<2; k++) { if ( !mirror_coordinates && ( ((int)x+(int)i < 0) || ((int)x+(int)i > input_cube->wX-1) || ((int)y+(int)j < 0) || ((int)y+(int)j > input_cube->wY-1) || ((int)z+(int)k < 0) || ((int)z+(int)k > input_cube->wZ-1)) ) {//coord out of input_cube ->ignore value //radius_values[ind_val] = (float)distMap->datac((unsigned short int)z, (unsigned short int)y, (unsigned short int)x) + 1; //min radius cannot be bigger than @ location in map + 1 //ind_val++; } else {//valid coord float distance_temp; distance_temp = (float)distMap->datac(\ coord_trans_int((int)floor(z+k), input_cube->wZ), \ coord_trans_int((int)floor(y+j), input_cube->wY), \ coord_trans_int((int)floor(x+i), input_cube->wX) ) \ * input_cube->get_spot(\ coord_trans_int((int)floor(z+k), input_cube->wZ), \ coord_trans_int((int)floor(y+j), input_cube->wY), \ coord_trans_int((int)floor(x+i), input_cube->wX) ); //radius = 0 if voxel is PORE float x_fracpart=0, y_fracpart=0, z_fracpart=0; // = fracpart if x>0, =1-fracpart if x<0 x_fracpart = x - (float)(int)floor(x); if (i>0) x_fracpart = (float)1 - x_fracpart; y_fracpart = y - (float)(int)floor(y); if (j>0) y_fracpart = (float)1 - y_fracpart; z_fracpart = z - (float)(int)floor(z); if (k>0) z_fracpart = (float)1 - z_fracpart; //if (i>0) // x_fracpart = (float)1 - x_fracpart; //else // x_fracpart = x - (float)(int)floor(x); //if (j>0) // y_fracpart = (float)1 - y_fracpart; //else // y_fracpart = y - (float)(int)floor(y); //if (k>0) // z_fracpart = (float)1 - z_fracpart; //else // z_fracpart = z - (float)(int)floor(z); //if (i>0) // x_fracpart = 1-(coord_trans_float(x, input_cube->wX)-(int)coord_trans_float(x, input_cube->wX)); //else // x_fracpart = (coord_trans_float(x, input_cube->wX)-(int)coord_trans_float(x, input_cube->wX)); //if (j>0) // y_fracpart = 1-(coord_trans_float(y, input_cube->wY)-(int)coord_trans_float(y, input_cube->wY)); //else // y_fracpart = (coord_trans_float(y, input_cube->wY)-(int)coord_trans_float(y, input_cube->wY)); //if (k>0) // z_fracpart = 1-(coord_trans_float(z, input_cube->wZ)-(int)coord_trans_float(z, input_cube->wZ)); //else // z_fracpart = coord_trans_float(z, input_cube->wZ)-(int)coord_trans_float(z, input_cube->wZ); if (distance_temp == 0) //simple Pythagore { //radius_values[ind_val] = sqrt( pow(x_fracpart,2)+pow(y_fracpart,2)+pow(z_fracpart,2) ); //distance_temp = sqrt( pow(x_fracpart,2)+pow(y_fracpart,2)+pow(z_fracpart,2) ); if ( (i==0) + (j==0) + (k==0) == 3) continue; //center voxel -> do nothing else distance_temp = sqrt( (i!=0)*pow(x_fracpart,2) \ + (j!=0)*pow(y_fracpart,2) \ + (k!=0)*pow(z_fracpart,2) ); } else switch ((i==0) + (j==0) + (k==0)) { case 2: if (i) //radius_values[ind_val] = x_fracpart; distance_temp += x_fracpart; else if (j) //radius_values[ind_val] = y_fracpart; distance_temp += y_fracpart; else if (k) //radius_values[ind_val] = z_fracpart; distance_temp += z_fracpart; //radius_values[ind_val] += distance_temp; break; case 1: if (i==0) { //determine coord closest to convolution voxel, then calc min possible radius (Pythagore) if ( y_fracpart < z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(y_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(y_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); else //z closest //radius_values[ind_val] = sqrt( pow(z_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(z_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); } else if (j==0) { if ( x_fracpart < z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(x_fracpart+1,2)+pow(z_fracpart+distance_temp-1,2) ); else //z closest //radius_values[ind_val] = sqrt( pow(z_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(z_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); } else if (k==0) { if ( x_fracpart < y_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(x_fracpart+1,2)+pow(y_fracpart+distance_temp-1,2) ); else //y closest //radius_values[ind_val] = sqrt( pow(y_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); distance_temp = sqrt( pow(y_fracpart+1,2)+pow(x_fracpart+distance_temp-1,2) ); } break; case 0: //determine coord closest to convolution voxel, then calc if (distance_temp == 1) { if ( x_fracpart <= y_fracpart ) //x closest { if ( x_fracpart <= z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart+1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); } else //y closest { if ( y_fracpart <= z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+1, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+1, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+1, 2) ); } } else { if ( x_fracpart <= y_fracpart ) //x closest { if ( x_fracpart <= z_fracpart ) //x closest //radius_values[ind_val] = sqrt( pow(x_fracpart+distance_temp-1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart+distance_temp-1, 2) + pow(y_fracpart, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); } else //y closest { if ( y_fracpart <= z_fracpart ) //y closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+distance_temp-1, 2) + pow(z_fracpart, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart+distance_temp-1, 2) + pow(z_fracpart, 2) ); else //z closest //radius_values[ind_val] = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); distance_temp = sqrt( pow(x_fracpart, 2) + pow(y_fracpart, 2) + pow(z_fracpart+distance_temp-1, 2) ); } } break; default: //case 3: //center voxel, ind_val= -> do nothing //radius_values[ind_val] = 9999999; //ignore center voxel by setting large value -> not min break; } //ind_val++; if (r > distance_temp) r = distance_temp; //get minimum value for r } } //r = mini(radius_values, 27); //get min of the min radius //added for debug: for empty samples, need to limit size of jump in FPT //r = min(r, max_jump_size*lambda_pix); or lambdaForFiber_pix ?? //delete [] radius_values; //if (r <= FPT_boundary_layer_lambdaForFiber_pix) //<= *=* to take care of rod=0 case -> still pool the exp distribution, not enough, cause lambdaForFiber_pix=0 in this case if (r <= FPT_boundary_layer_lambdaForFiber_pix) //fpt_flag = true; //while (r > FPT_boundary_layer_lambdaForFiber_pix) { //do r = -lambdaForFiber_pix*(float)log(1-random_number_gen.randExc()); //fpt_flag = false; //dist_n_disp_sq_fout << r << endl; //while (r >= 1); } else { fpt_flag = true; //r -= FPT_boundary_layer_lambdaForFiber_pix_4_fifth; //not to land at the interface - small likelyhood } } else { r = -lambdaForFiber_pix*(float)log( 1-random_number_gen.randExc() ); //r = min(r, max_jump_size*lambda_pix); or lambdaForFiber_pix ?? // to limit size of jump } //jumpinFiber ( ); {//tracer in FIBER float alpha, beta, gamma, mod; float dx, dy, dz; {//calc new tracer coord do //hypercube rejection method { alpha = 1-2*(float)random_number_gen.randExc();//randMT(1); beta = 1-2*(float)random_number_gen.randExc();//randMT(1); gamma = 1-2*(float)random_number_gen.randExc();//randMT(1); mod = pow(alpha, 2) + pow(beta, 2) + pow(gamma, 2); } while ( (mod < 0.001) || (mod > 1) ); mod = sqrt(mod); //switch (permut%3) //{ //case 0: dx = r*(alpha/mod); dy = r*(beta/mod); //jump coord dz = r*(gamma/mod); // break; //case 1: // dz = r*(alpha/mod); // dx = r*(beta/mod); //jump coord // dy = r*(gamma/mod); // break; //default: //case 2: // dy = r*(alpha/mod); // dz = r*(beta/mod); //jump coord // dx = r*(gamma/mod); // break; //} //permut ++; x2 = dx + x; y2 = dy + y; //hence new coordinates z2 = dz + z; } if (!mirror_coordinates) //check tracer didnt exit the roi { if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax) ) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) if ( (pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2)) > sq_MAX_ROI_RADIUS ) {//if tracer new coord out of sub vol (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub vol } } } //if ( (!fpt_flag) && (!input_cube->get_spot( if ( (!fpt_flag) && (!input_cube->get_spot( \ coord_trans_int((int)floor(z2), input_cube->wZ), \ coord_trans_int((int)floor(y2), input_cube->wY), \ coord_trans_int((int)floor(x2), input_cube->wX) ) ) ) //tracer (new coord) belongs to PORE space {//should only happen in *real/pure* random walk case hence the !fpt_flag #if RW_DO_DEBUG if ( fpt_flag ) // - to debug {//should not happen, not in FPT, unless wrong calc of r //may happen because of large jump with /sqrt_rod //may happen if alpha beta gamma randomly chosen such as 2 are <<1 and 1 is >>1 write_to_log("error FPT -3-"); } #endif //need to check if we skipped some PORE pixels in between, and locate the interface //check for every voxel on the path between origin (x)-fiber and end (x2)-pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final; float xi, yi, zi; //interface coord (should be int for at least one of them) { r_index_final = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x2 > x ) { for (x_index = (int)ceil(x); x_index <= (int)floor(x2); x_index++) //start @ x+1, finish @ x2 { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); // y = old + ratio * (new - old) z_index = (int)floor(z + r_index*(z2-z)); //check if coord valid in input_cube if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if pore r_index_final = r_index; break; } } } else { for (x_index = (int)floor(x); x_index >= (int)ceil(x2); x_index--) //start @ x, finish @ x2+1 { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); z_index = (int)floor(z + r_index*(z2-z)); if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final = r_index; break; } } } //for Y if ( y2 > y ) { for (y_index = (int)ceil(y); y_index <= (int)floor(y2); y_index++) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (y_index = (int)floor(y); y_index >= (int)ceil(y2); y_index--) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } //for Z if ( z2 > z ) { for (z_index = (int)ceil(z); z_index <= (int)floor(z2); z_index++) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (z_index = (int)floor(z); z_index >= (int)ceil(z2); z_index--) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } //if (r_index_final != (float)1) //re calc the coord to be at the interface - coord(i) xi = (x + r_index_final*(x2-x)); yi = (y + r_index_final*(y2-y)); zi = (z + r_index_final*(z2-z)); } //no reflection in this case, will always enter PORE space - no: reflection back to fiber possible if rod > 1 //original tracer coords (in FIBER space) before random jump float x3, y3, z3; //new coord after taking into account the interface coords float x4=0, y4=0, z4=0; //new coord after taking into account the possible jump over fiber if ( (float)random_number_gen.randExc() < inv_sqrt_rod ) { //fiber -> enter pore x3 = xi + inv_sqrt_rod * (x2-xi); y3 = yi + inv_sqrt_rod * (y2-yi); //going from LOW mobility media to HIGH mobility media (rod < 1 case, otherwise opposite) z3 = zi + inv_sqrt_rod * (z2-zi); //need to check if we skipped some fiber voxels //if (rod<1) { //need to check we didnt jump over some fiber voxels (rod<1 case - not necessarily can happen when rod>1 too) //check for every voxel on the path between origin (xi)-fiber and end (x3)-fiber or pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final2; float xi2, yi2, zi2; //interface coord (should be int for at least one of them) { r_index_final2 = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x3 > xi ) { for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) //start @ xi+1, finish @ x3 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) z_index = (int)floor(zi + r_index*(z3-zi)); //check if coord valid in input_cube if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if fiber r_index_final2 = r_index; break; } } } } else { for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) //start @ xi, finish @ x3+1 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final2 = r_index; break; } } } } //for Y if ( y3 > yi ) { for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } //for Z if ( z3 > zi ) { for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } if ( (r_index_final2 != (float)1) )//&& (r_index_final2 != (float)0) ) { //re calc the coord to be at the interface - coord(i) xi2 = (xi + r_index_final2*(x3-xi)); yi2 = (yi + r_index_final2*(y3-yi)); zi2 = (zi + r_index_final2*(z3-zi)); //then jump close to interface #2 minus a small distance (lambda_fiber?) r_index = sqrt( pow(xi2-xi,2) + pow(yi2-yi,2) + pow(zi2-zi,2) ); //distance between first interface and second one if (r_index > two_FPT_boundary_layer_lambdaForFiber_pix) //after jump should be at least lambda away from r_index -= FPT_boundary_layer_lambdaForFiber_pix;//FPT_boundary_layer_lambdaForFiber_pix_4_fifth; else r_index *= (float)0.55; //r_index = min(r_index, ); //r_index = min(r_index, ); x4 = xi2 + r_index * (xi - xi2); y4 = yi2 + r_index * (yi - yi2); z4 = zi2 + r_index * (zi - zi2); } else {//no jump over fiber, coord are valid, keep them x4 = x3; y4 = y3; z4 = z3; } } //else // rod>1 //{//no jump over fiber, coord are valid, keep them // x4 = x3; // y4 = y3; // z4 = z3; //} //re-calculate r because different than previously calculated (new=x4, previous before jump =x) r = sqrt( pow(x4-x, 2) + pow(y4-y, 2) + pow(z4-z, 2)); distance += PIXEL_SIZE*r; //in microns //save new coords as being current ones x = x4; y = y4; z = z4; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } } else { // reflected back (fiber -> fiber) //reflectFlag = 1; //{//bounce of how much would have entered freely // x3 = xi - (x2-xi); // y3 = yi - (y2-yi); //can cause problems (large bounce) // z3 = zi - (z2-zi); //} //{//bounce half way till the interface // x3 = x + (xi-x)/2; // y3 = y + (yi-y)/2; // z3 = z + (zi-z)/2; //} {//bounce back a random distance from the interface, depending on fractional_part(r) float fractpart, intpart; fractpart = modf (r , &intpart); //if = 0 or 1(0.999) then problem! can be pore (with x being integer) - very unlikely - taken care of with the 'continue" below x3 = x + (xi-x)*fractpart; y3 = y + (yi-y)*fractpart; z3 = z + (zi-z)*fractpart; } if ( !input_cube->get_spot( \ (unsigned short int)coord_trans_int((int)floor(z3), input_cube->wZ), \ (unsigned short int)coord_trans_int((int)floor(y3), input_cube->wY), \ (unsigned short int)coord_trans_int((int)floor(x3), input_cube->wX) ) ) {//re-pool the exp distribution #if RW_DO_DEBUG write_to_log("error reflection1 ->pore instead of fiber, will now 'continue'"); //debug #endif continue; } //re-calculate r because shorter than previously calculated (new=x3, previous before jump =x) //NOT taking into account the bouncing on the interface: //r = sqrt( pow(x3-x, 2)+pow(y3-y, 2)+pow(z3-z, 2) ); //taking into account the bouncing on the interface: //-> 150% distance between x and interface (sqrt(1.5)=1.224745) //r = (float)1.224745*sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ); //150%, bounce 1/2 way //taking into account the bouncing on the interface: r = sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ) \ + sqrt( pow(xi-x3, 2)+pow(yi-y3, 2)+pow(zi-z3, 2) ); // ||x,xi||+||xi,x3|| //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug //{//should not happen, not in FPT, unless wrong calc of r // #if RW_DO_DEBUG // write_to_log("error FPT -2-"); // #endif // //distance += pow(r, 2)/(2*lambda_pix); //in pixels // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns // distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //} //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns //new coord (x3) *should* be valid because x2 is (2nd reflection case, when bounce half way) //-> no, because of the reflection formula (the first one) //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // { // //write_to_log("simultaneous diffusion: pore->fiber -- ERROR2 -- "); // //if tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} //save new coords as being current ones x = x3; y = y3; z = z3; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } //{//need to check if we passed some FIBER pixels in between // //check for every voxel on the path // //between origin (xi) at interface and end (x3) voxels - stop as soon as reach FIBER // float r_index_final; // int x_index, y_index, z_index; // float r_index; // r_index = (float)1; // r_index_final = (float)1; // //for X // if ( x3 > xi ) // { // for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) //tracer belongs to FIBER space // {//if fiber // r_index_final = r_index; // break; // } // } // } // else // { // for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index-1, input_cube->wX) ) ) // { // r_index_final = r_index;//((float)(x_index+1) - xi) / (x3 - xi); // break; // } // } // } // //for Y // if ( y3 > yi ) // { // for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index); // break; // } // } // } // else // { // for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index-1, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index);//((float)(y_index+1) - yi) / (y3 - yi) ); // break; // } // } // } // //for Z // if ( z3 > zi ) // { // for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index); // break; // } // } // } // else // { // for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index-1, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index);//((float)(z_index+1) - zi) / (z3 - zi) ); // break; // } // } // } // if (r_index_final != (float)1) // {//re calc the coord to be at the interface // r_index_final -= (float)0.001; //to be still in pore space, just before the interface with fiber // x3 = (xi + r_index_final*(x3-xi)); // y3 = (yi + r_index_final*(y3-yi)); // z3 = (zi + r_index_final*(z3-zi)); // } // //else, (x3,y3,z3) is ok //} ////re-calculate r because LONGER than previously calculated (new=x3, previous before jump =x) //r = sqrt( pow(x3-x, 2) + pow(y3-y, 2) + pow(z3-z, 2)); //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // //if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // if ( pow(x3-xo, 2) + pow(y3-yo, 2) + pow(z3-zo, 2) > sq_MAX_ROI_RADIUS) // {//if tracer new coord out of sub cube (ROI) // //happens because of 1/sqrt_rod is a LARGE number // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} //{ // //calc distance // distance += PIXEL_SIZE*r; //in microns // //set new coords // x = x3; // y = y3; // z = z3; //} } else //tracer (new coord after jump) belongs to FIBER space { if ( fpt_flag ) {//keep new coordinates, tracer did NOT jump over some pore voxels x = x2; y = y2; z = z2; //no need to check if coords valid was done above - after //calc new tracer coord distance += PIXEL_SIZE_sq_over_2lambdaForFiber_micron * pow(r, 2); } else { //need to check tracer didnt jump over PORE voxels //check for every voxel on the path between origin (x)-fiber and end (x2)-pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final; float xi, yi, zi; //interface coord (should be int for at least one of them) { r_index_final = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x2 > x ) { for (x_index = (int)ceil(x); x_index <= (int)floor(x2); x_index++) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); // y = old + ratio * (new - old) z_index = (int)floor(z + r_index*(z2-z)); if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) //added y if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if pore r_index_final = r_index; break; } } } else { for (x_index = (int)floor(x); x_index >= (int)ceil(x2); x_index--) { r_index = ((float)x_index - x) / (x2 - x); y_index = (int)floor(y + r_index*(y2-y)); z_index = (int)floor(z + r_index*(z2-z)); if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final = r_index;//((float)(x_index+1) - x) / (x2 - x); break; } } } //for Y if ( y2 > y ) { for (y_index = (int)ceil(y); y_index <= (int)floor(y2); y_index++) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (y_index = (int)floor(y); y_index >= (int)ceil(y2); y_index--) { r_index = ((float)y_index - y) / (y2 - y); x_index = (int)floor(x + r_index*(x2-x)); z_index = (int)floor(z + r_index*(z2-z)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(y_index+1) - y) / (y2 - y) ); break; } } } //for Z if ( z2 > z ) { for (z_index = (int)ceil(z); z_index <= (int)floor(z2); z_index++) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index); break; } } } else { for (z_index = (int)floor(z); z_index >= (int)ceil(z2); z_index--) { r_index = ((float)z_index - z) / (z2 - z); x_index = (int)floor(x + r_index*(x2-x)); y_index = (int)floor(y + r_index*(y2-y)); if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) if (!input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final = min(r_index_final, r_index);//((float)(z_index+1) - z) / (z2 - z) ); break; } } } if (r_index_final != (float)1) {//means tracer jumped over PORE voxels //should only happen in *real/pure* random walk case #if RW_DO_DEBUG if ( fpt_flag ) // - to debug {//should not happen, not in FPT, unless wrong calc of r may happen because of large jump with /sqrt_rod write_to_log("error FPT -4-"); } #endif //re calc the coord to be at the interface - coord(i) xi = (x + r_index_final*(x2-x)); yi = (y + r_index_final*(y2-y)); zi = (z + r_index_final*(z2-z)); //no reflection in this case, will always enter PORE space - not true if rod>1 //original tracer coords (in FIBER space) before random jump float x3, y3, z3; //new coord after taking into account the interface coords float x4=0, y4=0, z4=0; //new coord after taking into account the possible jump over fiber if ( (float)random_number_gen.randExc() < inv_sqrt_rod ) { //fiber -> enter pore { x3 = xi + inv_sqrt_rod * (x2-xi); y3 = yi + inv_sqrt_rod * (y2-yi); //going from LOW mobility media to HIGH mobility media z3 = zi + inv_sqrt_rod * (z2-zi); } //need to check if we skipped some fiber voxels //if (rod<1) { //need to check we didnt jump over some fiber voxels (rod<1 case - not necessarily can happen when rod>1 too) //check for every voxel on the path between origin (xi)-fiber and end (x3)-fiber or pore voxels //- stop as soon as we meet a pore voxel, to get the corresponding interface coord float r_index_final2; float xi2, yi2, zi2; //interface coord (should be int for at least one of them) { r_index_final2 = (float)1; int x_index, y_index, z_index; float r_index; r_index = (float)1; //for X if ( x3 > xi ) { for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) //start @ xi+1, finish @ x3 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) z_index = (int)floor(zi + r_index*(z3-zi)); //check if coord valid in input_cube if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) {//if fiber r_index_final2 = r_index; break; } } } } else { for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) //start @ xi, finish @ x3+1 { if ((float)x_index != xi) { r_index = ((float)x_index - xi) / (x3 - xi); y_index = (int)floor(yi + r_index*(y3-yi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index-1, input_cube->wX) ) ) { r_index_final2 = r_index; break; } } } } //for Y if ( y3 > yi ) { for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) { if ((float)y_index != yi) { r_index = ((float)y_index - yi) / (y3 - yi); x_index = (int)floor(xi + r_index*(x3-xi)); z_index = (int)floor(zi + r_index*(z3-zi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (z_index >= (int)floor(min(zi, z3))) && (z_index < (int)ceil(max(zi, z3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index-1, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } //for Z if ( z3 > zi ) { for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } else { for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) { if ((float)z_index != zi) { r_index = ((float)z_index - zi) / (z3 - zi); x_index = (int)floor(xi + r_index*(x3-xi)); y_index = (int)floor(yi + r_index*(y3-yi)); if ((x_index >= (int)floor(min(xi, x3))) && (x_index < (int)ceil(max(xi, x3))) && (y_index >= (int)floor(min(yi, y3))) && (y_index < (int)ceil(max(yi, y3))) ) if (input_cube->get_spot( \ coord_trans_int(z_index-1, input_cube->wZ), \ coord_trans_int(y_index, input_cube->wY), \ coord_trans_int(x_index, input_cube->wX) ) ) { r_index_final2 = min(r_index_final2, r_index); break; } } } } if (r_index_final2 != (float)1) { //re calc the coord to be at the interface - coord(i) xi2 = (xi + r_index_final2*(x3-xi)); yi2 = (yi + r_index_final2*(y3-yi)); zi2 = (zi + r_index_final2*(z3-zi)); //then jump close to interface #2 minus a small distance (lambda_fiber?) r_index = sqrt( pow(xi2-xi,2) + pow(yi2-yi,2) + pow(zi2-zi,2) ); //distance between first interface and second one if (r_index > two_FPT_boundary_layer_lambdaForFiber_pix) //FPT_boundary_layer_lambdaForFiber_pix) //after jump should be at least lambda away from r_index -= FPT_boundary_layer_lambdaForFiber_pix;//FPT_boundary_layer_lambdaForFiber_pix_4_fifth; else r_index *= (float)0.55; //not good enough alone, ok with above statement //r_index = min(r_index, ); //r_index = min(r_index, ); x4 = xi2 + r_index * (xi - xi2); y4 = yi2 + r_index * (yi - yi2); z4 = zi2 + r_index * (zi - zi2); } else {//no jump over fiber, coord are valid, keep them x4 = x3; y4 = y3; z4 = z3; } } //else // rod>1 //{//no jump over fiber, coord are valid, keep them // x4 = x3; // y4 = y3; // z4 = z3; //} //re-calculate r because different than previously calculated (new=x4, previous before jump =x) r = sqrt( pow(x4-x, 2) + pow(y4-y, 2) + pow(z4-z, 2)); distance += PIXEL_SIZE*r; //in microns //save new coords as being current ones x = x4; y = y4; z = z4; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } } else { // reflected back (fiber -> fiber) //reflectFlag = 1; //{//bounce of how much would have entered freely // x3 = xi - (x2-xi); // y3 = yi - (y2-yi); //can cause problems (large bounce) // z3 = zi - (z2-zi); //} //{//bounce half way till the interface // x3 = x + (xi-x)/2; // y3 = y + (yi-y)/2; // z3 = z + (zi-z)/2; //} {//bounce back a random distance from the interface, depending on fractional_part(r) float fractpart, intpart; fractpart = modf (r , &intpart); //if = 0 or 1 then problem! could enter pore (with x being integer) - very unlikely - taken care of with the "continue" below x3 = x + (xi-x)*fractpart; y3 = y + (yi-y)*fractpart; z3 = z + (zi-z)*fractpart; } if ( !input_cube->get_spot( \ (unsigned short int)coord_trans_int((int)floor(z3), input_cube->wZ), \ (unsigned short int)coord_trans_int((int)floor(y3), input_cube->wY), \ (unsigned short int)coord_trans_int((int)floor(x3), input_cube->wX) ) ) {//re-pool the exp distribution #if RW_DO_DEBUG write_to_log("error reflection1 ->pore instead of fiber, will now 'continue'"); //debug #endif continue; } //re-calculate r because shorter than previously calculated (new=x3, previous before jump =x) //NOT taking into account the bouncing on the interface: //r = sqrt( pow(x3-x, 2)+pow(y3-y, 2)+pow(z3-z, 2) ); //taking into account the bouncing on the interface: //-> 150% distance between x and interface (sqrt(1.5)=1.224745) //r = (float)1.224745*sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ); //150%, bounce 1/2 way //taking into account the bouncing on the interface: r = sqrt( pow(xi-x, 2)+pow(yi-y, 2)+pow(zi-z, 2) ) \ + sqrt( pow(xi-x3, 2)+pow(yi-y3, 2)+pow(zi-z3, 2) ); // ||x,xi||+||xi,x3|| //calc distance //if ( fpt_flag ) //(r>lambda_pix) && fpt_flag ) - to debug //{//should not happen, not in FPT, unless wrong calc of r // #if RW_DO_DEBUG // write_to_log("error FPT -2-"); // #endif // //distance += pow(r, 2)/(2*lambda_pix); //in pixels // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambda_micron); //in microns // distance += PIXEL_SIZE_sq_over_2lambda_micron * pow(r, 2); //in microns //} //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns //new coord (x3) *should* be valid because x2 is (2nd reflection case, when bounce half way) //-> no, because of the reflection formula (the first one) //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // { // //write_to_log("simultaneous diffusion: pore->fiber -- ERROR2 -- "); // //if tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} //save new coords as being current ones x = x3; y = y3; z = z3; if (!mirror_coordinates) {//check coords are valid if (!sub_ROI_is_sphere) {//this is a cube type of sub-volume if ((x2 < xROImin) || (x2 >= xROImax) || (y2 < yROImin) || (y2 >= yROImax) || (z2 < zROImin) || (z2 >= zROImax )) break; //tracer exited the ROI sub cube } else {//this is a sphere type of sub-volume //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // > sq_MAX_ROI_RADIUS) //if (!random_starting_point) { if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) //tracer new coord out of sub cube (ROI) //xf = x; //yf = y; //final coordinates just before exiting the ROI //zf = z; break; //tracer exited the ROI sub cube } //else //{//criteria to consider tracer has traveled enough // random starting point case /*if (abs(x2-xo) > ) break; if (abs(y2-yo) > ) break; if (abs(z2-zo) > ) break; */ //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) //break; //} } } } //{//need to check if we passed some FIBER pixels in between // //check for every voxel on the path // //between origin (xi) at interface and end (x3) voxels - stop as soon as reach FIBER // float r_index_final; // int x_index, y_index, z_index; // float r_index; // r_index = (float)1; // r_index_final = (float)1; // //for X // if ( x3 > xi ) // { // for (x_index = (int)ceil(xi); x_index <= (int)floor(x3); x_index++) // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // y = old + ratio * (new - old) // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) //added y // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) //tracer belongs to FIBER space // { // r_index_final = r_index; // break; // } // } // } // else // { // for (x_index = (int)floor(xi); x_index >= (int)ceil(x3); x_index--) // { // r_index = ((float)x_index - xi) / (x3 - xi); // y_index = (int)floor(yi + r_index*(y3-yi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) //added y // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index-1, input_cube->wX) ) ) // { // r_index_final = r_index;//((float)(x_index+1) - xi) / (x3 - xi); // break; // } // } // } // //for Y // if ( y3 > yi ) // { // for (y_index = (int)ceil(yi); y_index <= (int)floor(y3); y_index++) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index); // break; // } // } // } // else // { // for (y_index = (int)floor(yi); y_index >= (int)ceil(y3); y_index--) // { // r_index = ((float)y_index - yi) / (y3 - yi); // x_index = (int)floor(xi + r_index*(x3-xi)); // z_index = (int)floor(zi + r_index*(z3-zi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (z_index >= (int)floor(min(z, z2))) && (z_index < (int)ceil(max(z, z2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index-1, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index);//((float)(y_index+1) - yi) / (y3 - yi) ); // break; // } // } // } // //for Z // if ( z3 > zi ) // { // for (z_index = (int)ceil(zi); z_index <= (int)floor(z3); z_index++) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index); // break; // } // } // } // else // { // for (z_index = (int)floor(zi); z_index >= (int)ceil(z3); z_index--) // { // r_index = ((float)z_index - zi) / (z3 - zi); // x_index = (int)floor(xi + r_index*(x3-xi)); // y_index = (int)floor(yi + r_index*(y3-yi)); // if ((x_index >= (int)floor(min(x, x2))) && (x_index < (int)ceil(max(x, x2))) && // (y_index >= (int)floor(min(y, y2))) && (y_index < (int)ceil(max(y, y2))) ) // if (input_cube->get_spot( \ // coord_trans_int(z_index-1, input_cube->wZ), \ // coord_trans_int(y_index, input_cube->wY), \ // coord_trans_int(x_index, input_cube->wX) ) ) // { // r_index_final = min(r_index_final, r_index);//((float)(z_index+1) - zi) / (z3 - zi) ); // break; // } // } // } // if (r_index_final != (float)1) // {//re calc the coord to be at the interface // r_index_final -= (float)0.001; //to be still in pore space, just before the fiber space // x3 = (xi + r_index_final*(x3-xi)); // y3 = (yi + r_index_final*(y3-yi)); // z3 = (zi + r_index_final*(z3-zi)); // } // //else, (x3,y3,z3) is ok //} ////re-calculate r because LONGER than previously calculated (new=x3, previous before jump =x) //r = sqrt( pow(x3-x, 2) + pow(y3-y, 2) + pow(z3-z, 2)); //if (!sub_ROI_is_sphere) //{ // if ((x3 < xROImin) || (x3 >= xROImax) || // (y3 < yROImin) || (y3 >= yROImax) || // (z3 < zROImin) || (z3 >= zROImax) ) //this is a cube type of sub-volume // break; //tracer exited the ROI sub cube //} //else //{ // //if ( pow(x3-(float)x_start, 2) + pow(y3-(float)y_start, 2) + pow(z3-(float)z_start, 2) \ // // > sq_MAX_ROI_RADIUS) //this is a sphere type of sub-volume // if ( pow(x3-xo, 2) + pow(y3-yo, 2) + pow(z3-zo, 2) > sq_MAX_ROI_RADIUS) // {//if tracer new coord out of sub cube (ROI) // //happens because of 1/sqrt_rod is a LARGE number // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } //} //{ // //calc distance // distance += PIXEL_SIZE*r; //in microns // //set new coords // x = x3; // y = y3; // z = z3; // if (!mirror_coordinates) // {//check coords are valid // if (!sub_ROI_is_sphere) // {//this is a cube type of sub-volume // if ((x2 < xROImin) || (x2 >= xROImax) || // (y2 < yROImin) || (y2 >= yROImax) || // (z2 < zROImin) || (z2 >= zROImax )) // break; //tracer exited the ROI sub cube // } // else // {//this is a sphere type of sub-volume // //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // // > sq_MAX_ROI_RADIUS) // //if (!random_starting_point) // { // if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) // //tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } // //else // //{//criteria to consider tracer has traveled enough // random starting point case // /*if (abs(x2-xo) > ) // break; // if (abs(y2-yo) > ) // break; // if (abs(z2-zo) > ) // break; // */ // //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) // //break; // //} // } // } //} } else //didnt jump over PORE {//keep new coordinates x = x2; y = y2; z = z2; //no need to check coords was done above //if (!mirror_coordinates) //{//check coords are valid // if (!sub_ROI_is_sphere) // {//this is a cube type of sub-volume // if ((x2 < xROImin) || (x2 >= xROImax) || // (y2 < yROImin) || (y2 >= yROImax) || // (z2 < zROImin) || (z2 >= zROImax )) // break; //tracer exited the ROI sub cube // } // else // {//this is a sphere type of sub-volume // //if ( pow(x2-(float)x_start, 2) + pow(y2-(float)y_start, 2) + pow(z2-(float)z_start, 2) \ // // > sq_MAX_ROI_RADIUS) // //if (!random_starting_point) // { // if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > sq_MAX_ROI_RADIUS) // //tracer new coord out of sub cube (ROI) // //xf = x; // //yf = y; //final coordinates just before exiting the ROI // //zf = z; // break; //tracer exited the ROI sub cube // } // //else // //{//criteria to consider tracer has traveled enough // random starting point case // /*if (abs(x2-xo) > ) // break; // if (abs(y2-yo) > ) // break; // if (abs(z2-zo) > ) // break; // */ // //if ( pow(x2-xo, 2) + pow(y2-yo, 2) + pow(z2-zo, 2) > max_displacement_diag_pixels_sq) // //break; // //} // } //} //calc distance //if ( fpt_flag ) //(r>lambdaForFiber_pix) && (!dlg->CHECK_TRUE_RANDOM_WALK.GetCheck()) ) // //distance += pow(r, 2)/(2*lambdaForFiber); // //distance += PIXEL_SIZE_sq*pow(r, 2)/(2*lambdaForFiber_micron); // distance += PIXEL_SIZE_sq_over_2lambdaForFiber_micron * pow(r, 2); //else //distance += r; //in pixels distance += PIXEL_SIZE*r; //in microns } } } } } if (fpt_flag) ++num_jumps_fiber_fpt; else ++num_jumps_fiber_nonfpt; } //simulationResults (); if (distance >= dlg->field_rw_MAX_DISTANCE_float) { //write_to_log("max_distance exceeded"); max_distance_exceeded = true; max_distance_exceeded_count ++; //dist_n_disp_sq_fout << "max_distance exceeded" << endl; //not to write it for each and every tracer continue; //go to next tracer } else { if (distance >= dataCollectionPoint) { //float dispx, dispy, dispz; distance_structure temp_entry; //dataCollectionPoint = distance + (float)dataCollectionPoint_increment; //fmod dataCollectionPoint = dataCollectionPoint_increment * (float)(int)( 1 + distance / dataCollectionPoint_increment); count_distances ++; temp_entry.distance = distance; //in microns already //temp_entry.dispx = pow(x-(float)xo, 2); //temp_entry.dispy = pow(y-(float)yo, 2); //in pixels //temp_entry.dispz = pow(z-(float)zo, 2); temp_entry.dispx = PIXEL_SIZE_sq*pow(x-xo, 2); temp_entry.dispy = PIXEL_SIZE_sq*pow(y-yo, 2); //in microns temp_entry.dispz = PIXEL_SIZE_sq*pow(z-zo, 2); if ( dlg->CHECK_SAVE_DETAILED_TRACE.GetCheck() ) {//save results to file (debug?) dist_n_disp_sq_fout << temp_entry.distance << ","; dist_n_disp_sq_fout << temp_entry.dispx << ","; dist_n_disp_sq_fout << temp_entry.dispy << ","; dist_n_disp_sq_fout << temp_entry.dispz << ","; //for debug dist_n_disp_sq_fout << x << "," << y << "," << z << "," ; //dist_n_disp_sq_fout << input_cube->get_spot( (unsigned short int)z, (unsigned short int)y, (unsigned short int)x ) << endl;//<< "," ; //dist_n_disp_sq_fout << (int)distMap->datac((unsigned short int)z, (unsigned short int)y, (unsigned short int)x) << endl; dist_n_disp_sq_fout << endl; //save tracer position into path map ds_cube_tracer_path->set_spot_1( (unsigned short int)coord_trans_float(z, input_cube->wZ), (unsigned short int)coord_trans_float(y, input_cube->wY), (unsigned short int)coord_trans_float(x, input_cube->wX) ); } distance_structure_list_ptr->insertAtEnd(temp_entry); } } //if (distance >= max_allowed_distance_pix) //{//consider tracer lost inside closed pore -> should not happen often in simultaneous diffusion case // //erase all previous values in the distance_structure_list_ptr // particle --; // continue; //} //update progress bar //if (particle%2 == 0) dlg->pbar_simulation_progress.SetPos((short int)((float)particle/(float)total_particle_per_roi*(float)total_particle_for_progressbar)); //dlg->pbar_simulation_progress.StepIt(); //dist_n_disp_sq_fout << endl; //to separate results for each particle if (distance > max_distance) max_distance = distance; ////save tracer position into path map //ds_cube_tracer_path->set_spot_1( (unsigned short int)z, (unsigned short int)y, (unsigned short int)x ); } if ( dlg->CHECK_SAVE_DETAILED_TRACE.GetCheck() ) { // Write tracer path result to file char tracer_path_file_out [1024]; sprintf(tracer_path_file_out, "%s_tracer#%d_path", init_file_name_for_radius, particle); ds_cube_tracer_path->writeToFile( tracer_path_file_out, true ); delete ds_cube_tracer_path; } int intermediate_cv_check; intermediate_cv_check = min( 200, (int)((float)cv_num_tracers/(float)5) ); if ( check_cv & ( (particle%cv_num_tracers == 0) | (particle%intermediate_cv_check == 0) ) ) {//test convergence criteria here if ( (rod!=0) | (total_file_number==1) ) //ie rod!= 0 OR rod=0 and 1 single sub cube - mirroring case //ok to proceed {//do stats if (particle%cv_num_tracers == 0) { write_to_log("**DATE**"); sprintf(win_txt, "check CV after # tracers, %d", particle); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); } // - least square fit - linear regression //http://mathworld.wolfram.com/LeastSquaresFitting.html //STEP1 //get bucket width: we want the buckets to cover the distance //betweeen 0 and a limited_max_distance, such that percent_max_distance of the datapoints are < limited_max_distance (and > 0 of course) //get limited_max_distance first: use a total_num_buckets buckets histogram to determine it // Grab start of list pointer LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); unsigned int total_num_buckets, bucket_index; //unsigned long long count_number_datapoints_in_distance_structure = 0; total_num_buckets = 20; //dlg->field_rw_NUM_BUCKETS1_int; //10000; //20 is enough unsigned long long total_count_distance_bucket = 0, total_count_distance_bucket2 = 0;; float percent_max_distance; percent_max_distance = dlg->field_rw_PERCENT_MAX_DISTANCE_float;//(float).5; float percent_densest_bucket; percent_densest_bucket = dlg->field_rw_PERCENT_DENSEST_BUCKET_float;//(float).9; float limited_max_distance; float bucket_width; bucket_width = (max_distance+(float)1)/(float)total_num_buckets; //in microns unsigned long long * count_distance_bucket1 = new unsigned long long [total_num_buckets+1]; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket1[bucket_index] = (unsigned long long)0; } while ( distance_structure_list_node_ptr ) { distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket1[bucket_index] ++; distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; //count_number_datapoints_in_distance_structure++; //for debug only } //write_to_log("# of datapoints %d", count_number_datapoints_in_distance_structure); //for debug only total_count_distance_bucket = 0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { total_count_distance_bucket += count_distance_bucket1[bucket_index]; } total_count_distance_bucket2 = 0; bucket_index=0; while (total_count_distance_bucket2 < percent_max_distance*total_count_distance_bucket) { total_count_distance_bucket2 += count_distance_bucket1[bucket_index]; bucket_index ++; } delete [] count_distance_bucket1; limited_max_distance = ((float)bucket_index+1)*bucket_width; //STEP2 total_num_buckets = dlg->field_rw_NUM_BUCKETS2_int;//20; double * distance_bucket = new double [total_num_buckets+1]; double * disp2x_bucket = new double [total_num_buckets+1]; double * disp2y_bucket = new double [total_num_buckets+1]; double * disp2z_bucket = new double [total_num_buckets+1]; double * sd_distance_bucket = new double [total_num_buckets+1]; double * sd_dispx_bucket = new double [total_num_buckets+1]; double * sd_dispy_bucket = new double [total_num_buckets+1]; double * sd_dispz_bucket = new double [total_num_buckets+1]; unsigned long long * count_distance_bucket = new unsigned long long [total_num_buckets+1]; unsigned short int num_used_buckets; double SSxx_distance, SSyy_dispx, SSyy_dispy, SSyy_dispz; //(SSxx and SSyys) double SSxy_ddx, SSxy_ddy, SSxy_ddz; //(SSxys) float av_distance, av_dispx, av_dispy, av_dispz; do { bucket_width = limited_max_distance/(float)total_num_buckets; float num_times_dataCollectionPoint_increment; num_times_dataCollectionPoint_increment = (float)1; if (bucket_width < num_times_dataCollectionPoint_increment*dataCollectionPoint_increment) { write_to_log("bucket_width @ STEP2 too small = %f", bucket_width); bucket_width = num_times_dataCollectionPoint_increment*dataCollectionPoint_increment; //to ensure enough write_to_log("increasing size to %f *dataCollectionPoint_increment = %f", num_times_dataCollectionPoint_increment, bucket_width); } unsigned long long max_num_in_a_bucket = 0; //unsigned long long counter; //count_distances //total_num_buckets = (unsigned short int)(limited_max_distance/bucket_width) + 1; //float av_buckets; //unsigned short int num_av_buckets; //av_buckets = (float)0; //num_av_buckets = 0; SSxx_distance = (double)0; SSyy_dispx = (double)0; SSyy_dispy = (double)0; SSyy_dispz = (double)0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket[bucket_index] = (unsigned long long)0; distance_bucket[bucket_index] = (double)0; disp2x_bucket[bucket_index] = (double)0; disp2y_bucket[bucket_index] = (double)0; disp2z_bucket[bucket_index] = (double)0; sd_distance_bucket[bucket_index] = (double)0; sd_dispx_bucket[bucket_index] = (double)0; sd_dispy_bucket[bucket_index] = (double)0; sd_dispz_bucket[bucket_index] = (double)0; } //LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket[bucket_index] ++; distance_bucket[bucket_index] += (double)temp_entry.distance; disp2x_bucket[bucket_index] += (double)temp_entry.dispx; disp2y_bucket[bucket_index] += (double)temp_entry.dispy; disp2z_bucket[bucket_index] += (double)temp_entry.dispz; } distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } // should depend on total_num_buckets? // depending on a fixed value -> not good //min_required_num_per_bucket = (unsigned int)(count_distances*0.02); //xx% of total # of datapoints -> not good either // depending on value in most populated bucket for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) if ( max_num_in_a_bucket < count_distance_bucket[bucket_index] ) max_num_in_a_bucket = count_distance_bucket[bucket_index]; min_required_num_per_bucket = (unsigned int)(max_num_in_a_bucket*percent_densest_bucket); //percent_densest_bucket% of most populated bucket num_used_buckets = 0; total_count_distance_bucket = 0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//summing everything - SS is just the overall SUM here -> temp variables to calc overall AVERAGES first if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) {//keep value only if there are enough datapoints in that particular bucket num_used_buckets ++; SSxx_distance += distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; total_count_distance_bucket += count_distance_bucket[bucket_index]; //SSxx_distance += distance_bucket[bucket_index] = ((double)bucket_index+.5) * (double)bucket_width; //SSxx_distance calc ???should??? involve middle of bucket value, instead of real data value SSyy_dispx += disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispy += disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispz += disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; } else { distance_bucket[bucket_index] = 0; disp2x_bucket[bucket_index] = 0; //averages per bucket disp2y_bucket[bucket_index] = 0; disp2z_bucket[bucket_index] = 0; } } percent_densest_bucket -= (float).01; //decrease % to ensure at least 3 buckets } while ( (num_used_buckets < 3) && (percent_densest_bucket >= 0) ); if (num_used_buckets == 1 ) { write_to_log("only 1 bucket!"); } //calc overall averages { av_distance = (float)SSxx_distance/(float)num_used_buckets; //xbar av_dispx = (float)SSyy_dispx/(float)num_used_buckets; av_dispy = (float)SSyy_dispy/(float)num_used_buckets; //ybar av_dispz = (float)SSyy_dispz/(float)num_used_buckets; } //now that we have the averages //calc standard errors for everything // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) { sd_distance_bucket[bucket_index] += pow((double)temp_entry.distance - distance_bucket[bucket_index], 2); sd_dispx_bucket[bucket_index] += pow((double)temp_entry.dispx - disp2x_bucket[bucket_index], 2); sd_dispy_bucket[bucket_index] += pow((double)temp_entry.dispy - disp2y_bucket[bucket_index], 2); sd_dispz_bucket[bucket_index] += pow((double)temp_entry.dispz - disp2z_bucket[bucket_index], 2); } } distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } //init variables { SSxx_distance = (double)0; //SSxx SSyy_dispx = (double)0; SSyy_dispy = (double)0; //SSyy SSyy_dispz = (double)0; SSxy_ddx = (double)0; SSxy_ddy = (double)0; //SSxy SSxy_ddz = (double)0; } //calc SS for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) if (distance_bucket[bucket_index] > 0) { SSxx_distance += pow(distance_bucket[bucket_index], 2); //SSxx=SUM(x^2) SSyy_dispx += pow(disp2x_bucket[bucket_index], 2); SSyy_dispy += pow(disp2y_bucket[bucket_index], 2); //SSyy=SUM(y^2) SSyy_dispz += pow(disp2z_bucket[bucket_index], 2); SSxy_ddx += distance_bucket[bucket_index] * disp2x_bucket[bucket_index]; SSxy_ddy += distance_bucket[bucket_index] * disp2y_bucket[bucket_index]; //SSxy=SUM(x*y) SSxy_ddz += distance_bucket[bucket_index] * disp2z_bucket[bucket_index]; sd_distance_bucket[bucket_index] = sqrt(sd_distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispx_bucket[bucket_index] = sqrt(sd_dispx_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispy_bucket[bucket_index] = sqrt(sd_dispy_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispz_bucket[bucket_index] = sqrt(sd_dispz_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); //write values to detailed cv data file rw_cv_data_detailed_fout << ", " << distance_bucket[bucket_index] << ", " \ << disp2x_bucket[bucket_index] << ", " \ << disp2y_bucket[bucket_index] << ", " \ << disp2z_bucket[bucket_index] << ", " \ << count_distance_bucket[bucket_index] << ", " << endl; } //delete temporary variables { delete [] count_distance_bucket; delete [] distance_bucket; delete [] disp2x_bucket; delete [] disp2y_bucket; delete [] disp2z_bucket; delete [] sd_distance_bucket; delete [] sd_dispx_bucket; delete [] sd_dispy_bucket; delete [] sd_dispz_bucket; } //save results for 0-intercept float Sumxx, Sumxyx, Sumxyy, Sumxyz; Sumxx = (float)SSxx_distance; Sumxyx = (float)SSxy_ddx; Sumxyy = (float)SSxy_ddy; Sumxyz = (float)SSxy_ddz; //finally: SSxx_distance = SSxx_distance - (double)num_used_buckets * pow((double)av_distance, 2); //SSxx=SUM(x^2)-n*xbar^2 SSyy_dispx = SSyy_dispx - (double)num_used_buckets * pow((double)av_dispx, 2); SSyy_dispy = SSyy_dispy - (double)num_used_buckets * pow((double)av_dispy, 2); //SSyy=SUM(y^2)-n*ybar^2 SSyy_dispz = SSyy_dispz - (double)num_used_buckets * pow((double)av_dispz, 2); SSxy_ddx = SSxy_ddx - (double)num_used_buckets * (double)av_distance * (double)av_dispx; SSxy_ddy = SSxy_ddy - (double)num_used_buckets * (double)av_distance * (double)av_dispy; //SSxy=SUM(x*y)-n*xbar*ybar SSxy_ddz = SSxy_ddz - (double)num_used_buckets * (double)av_distance * (double)av_dispz; //get a, b (y=a+bx) -- disp^2 vs distance //b=SSxy/SSxx float bx, by, bz; bx = (float)(SSxy_ddx/SSxx_distance); by = (float)(SSxy_ddy/SSxx_distance); bz = (float)(SSxy_ddz/SSxx_distance); //(all we want is the slope b) //now a=ybar-b*xbar float ax, ay, az; ax = av_dispx - bx*av_distance; ay = av_dispy - by*av_distance; az = av_dispz - bz*av_distance; float sigmax_distance, sigmay_dispx, sigmay_dispy, sigmay_dispz, covxy_ddx, covxy_ddy, covxy_ddz; //sigmax^2=SSxx/n sigmax_distance = (float)(SSxx_distance / num_used_buckets); //sigmay^2=SSyy/n sigmay_dispx = (float)(SSyy_dispx / (double)num_used_buckets); sigmay_dispy = (float)(SSyy_dispy / (double)num_used_buckets); sigmay_dispz = (float)(SSyy_dispz / (double)num_used_buckets); //Cov(x,y)=SSxy/n covxy_ddx = (float)(SSxy_ddx / (double)num_used_buckets); covxy_ddy = (float)(SSxy_ddy / (double)num_used_buckets); covxy_ddz = (float)(SSxy_ddz / (double)num_used_buckets); //s=sqrt( (SSyy-b*SSxy)/(n-2) ) float s_x, s_y, s_z; s_x = (float)sqrt( (SSyy_dispx-bx*SSxy_ddx)/((float)num_used_buckets-2) ); s_y = (float)sqrt( (SSyy_dispy-by*SSxy_ddy)/((float)num_used_buckets-2) ); s_z = (float)sqrt( (SSyy_dispz-bz*SSxy_ddz)/((float)num_used_buckets-2) ); //se(a)=s*sqrt( 1/n + xbar^2/SSxx ) float se_ax, se_ay, se_az; se_ax = s_x*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_ay = s_y*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_az = s_z*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); //se(b)=s/sqrt(SSxx) float se_bx, se_by, se_bz; se_bx = s_x / (float)sqrt(SSxx_distance); se_by = s_y / (float)sqrt(SSxx_distance); se_bz = s_z / (float)sqrt(SSxx_distance); //R^2=SSxy^2/SSxx/SSyy float R2_x, R2_y, R2_z; { R2_x = (float)(SSxy_ddx*SSxy_ddx/SSxx_distance/SSyy_dispx); R2_y = (float)(SSxy_ddy*SSxy_ddy/SSxx_distance/SSyy_dispy); R2_z = (float)(SSxy_ddz*SSxy_ddz/SSxx_distance/SSyy_dispz); } //now get D/D0 if (rod!=0) { D_D0x = (float)1.5*bx/lambda_micron; D_D0y = (float)1.5*by/lambda_micron; D_D0z = (float)1.5*bz/lambda_micron; } else {//rod == 0 D_D0x = (float)1.5*bx*ROI_porosity/lambda_micron; D_D0y = (float)1.5*by*ROI_porosity/lambda_micron; D_D0z = (float)1.5*bz*ROI_porosity/lambda_micron; } //stats_results_fout << "with zero (0) intercept" << endl; //stats_results_fout << ", " << 1.5*(float)(Sumxyx/Sumxx)/lambda_micron \ // << ", " << 1.5*(float)(Sumxyy/Sumxx)/lambda_micron \ // << ", " << 1.5*(float)(Sumxyz/Sumxx)/lambda_micron << endl; //now compare to previous values if (particle%cv_num_tracers == 0) { if ((unsigned long int)particle == cv_num_tracers) {//this is the first time D/D0 are calculated, can t compare them. store results for now //write_to_log("1st time D/D0 calculated"); D_D0x_previous = D_D0x; D_D0y_previous = D_D0y; D_D0z_previous = D_D0z; //write_to_log("D_D0x %f, D_D0y %f, D_D0z %f", D_D0x, D_D0y, D_D0z); //write_to_log("D/D0new D/D0old cv_criteria_sq < %f", cv_criteria_sq); //write_to_log("D/D0new D/D0old cv_criteria < %f", cv_criteria); rw_cv_data_fout << particle <<","<< D_D0x <<","<< D_D0y <<","<< D_D0z<<","; rw_cv_data_fout << "1,1,1," ; rw_cv_data_fout << R2_x <<","<< R2_y <<","<< R2_z << endl; } else {//now we can compare values to old ones, need ALL 3 (x,y,z) to match the criteria //write_to_log("D/D0new D/D0old cv_criteria_met"); float cv_temp_x, cv_temp_y, cv_temp_z; cv_temp_x = abs( (float)(D_D0x_previous-D_D0x) / (float)D_D0x_previous); cv_temp_y = abs( (float)(D_D0y_previous-D_D0y) / (float)D_D0y_previous); cv_temp_z = abs( (float)(D_D0z_previous-D_D0z) / (float)D_D0z_previous); //write_to_log("x, %f, %f, %f", D_D0x, D_D0x_previous, cv_temp_x); //write_to_log("y, %f, %f, %f", D_D0y, D_D0y_previous, cv_temp_y); //write_to_log("z, %f, %f, %f", D_D0z, D_D0z_previous, cv_temp_z); rw_cv_data_fout << particle <<","<< D_D0x <<","<< D_D0y <<","<< D_D0z <<","; rw_cv_data_fout << cv_temp_x <<","<< cv_temp_y <<","<< cv_temp_z <<","; rw_cv_data_fout << R2_x <<","<< R2_y <<","<< R2_z << endl; rw_converged = ( cv_temp_x < cv_criteria ); rw_converged &= ( cv_temp_y < cv_criteria ); rw_converged &= ( cv_temp_z < cv_criteria ); if (!rw_converged) { //save old values as new D_D0x_previous = D_D0x; D_D0y_previous = D_D0y; D_D0z_previous = D_D0z; //write_to_log("D_D0x, %f, D_D0y, %f, D_D0z, %f", D_D0x, D_D0y, D_D0z); } else { //rw_converged = true; //write_to_log("cved @ %f", cv_criteria); { total_particle_end_of_run += particle; //save total number of particles used rw_cv_data_fout << "cved @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_fout.close(); rw_cv_data_detailed_fout << particle <<","<< D_D0x <<","<< D_D0y <<","<< D_D0z<<","; rw_cv_data_detailed_fout << ",,," ; rw_cv_data_detailed_fout << R2_x <<","<< R2_y <<","<< R2_z << endl; rw_cv_data_detailed_fout << endl; rw_cv_data_detailed_fout << "cved @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_detailed_fout.close(); } } } //next check cv with at least cv_percent_new_data % more data cv_num_tracers = (unsigned long int)max(cv_num_tracers, (unsigned long int)((float)particle * cv_percent_new_data)); } else {//just write down data rw_cv_data_fout << particle <<","<< D_D0x <<","<< D_D0y <<","<< D_D0z<<","; rw_cv_data_fout << ",,," ; rw_cv_data_fout << R2_x <<","<< R2_y <<","<< R2_z << endl; } if (!rw_converged) {//detailed cv data rw_cv_data_detailed_fout << particle <<","<< D_D0x <<","<< D_D0y <<","<< D_D0z<<","; rw_cv_data_detailed_fout << ",,," ; rw_cv_data_detailed_fout << R2_x <<","<< R2_y <<","<< R2_z << endl; rw_cv_data_detailed_fout << endl; } } } particle++; //distance_structure_list_node_ptr = distance_structure_list_ptr->getTailPtr(); } while ( ( particle <= (int)total_particle_per_roi ) & !rw_converged & !dlg->stop_request); total_particle_end_of_run += particle-1; //save total number of particles used if (check_cv & !rw_converged) { //ie dlg->stop_request true or max # particles reached rw_cv_data_fout << "stopped @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_fout.close(); rw_cv_data_detailed_fout << "stopped @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_detailed_fout.close(); } if (max_distance_exceeded) { write_to_log("max_distance exceeded, check detailed trace coord. file"); max_distance_exceeded = false; } if (rod) //ie != 0 { //write down porosity from where tracers landed dist_n_disp_sq_fout << "ROI_porosity, " << ROI_porosity << endl; dist_n_disp_sq_fout << "porosity from tracer start, " << (float)count_tracer_start_hit_pore/(float)count_tracer_start_hit << endl; dist_n_disp_sq_fout << "num_jumps_pore_fpt, " << num_jumps_pore_fpt << ", "; dist_n_disp_sq_fout << "num_jumps_pore_nonfpt, " << num_jumps_pore_nonfpt << endl; dist_n_disp_sq_fout << "num_jumps_fiber_fpt, " << num_jumps_fiber_fpt << ", "; dist_n_disp_sq_fout << "num_jumps_fiber_nonfpt, " << num_jumps_fiber_nonfpt << endl; dist_n_disp_sq_fout << "total_particle_end_of_run, " << total_particle_end_of_run << ", "; dist_n_disp_sq_fout << "max_distance_exceeded_count, " << max_distance_exceeded_count; dist_n_disp_sq_fout << ",%, " << (float)max_distance_exceeded_count/(float)total_particle_end_of_run << endl; dist_n_disp_sq_fout << ",,displ_diag_exceeded_count, " << total_particle_end_of_run-max_distance_exceeded_count; dist_n_disp_sq_fout << ",%, " << (float)(total_particle_end_of_run-max_distance_exceeded_count)/(float)total_particle_end_of_run << endl; //close distance&displacement values file dist_n_disp_sq_fout << endl << "****" << endl << endl; dist_n_disp_sq_fout.close(); sprintf(win_txt, "Simulation: rep %d/%d (%d particles) done", file_number, total_file_number, particle-1); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); } else //rod==0 {//stats if (dlg->stop_request) { sprintf(win_txt, "stop_request received, attempting to do stats and finish now"); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); stats_results_per_file_number_fout << "CANCELED @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rod_zero_aborted_enough_data_for_stats = false; //rw_cv_data_fout.close(); //goto delvar; //goto dostats; } else { sprintf(win_txt, "Simulation: rep %d/%d (%d particles) done - now doing stats", file_number, total_file_number, total_particle_per_roi); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); } stats_results_per_file_number_fout << "num_jumps_pore_fpt, " << num_jumps_pore_fpt << ", "; stats_results_per_file_number_fout << "num_jumps_pore_nonfpt, " << num_jumps_pore_nonfpt << endl; stats_results_per_file_number_fout << "num_jumps_fiber_fpt, " << num_jumps_fiber_fpt << ", "; stats_results_per_file_number_fout << "num_jumps_fiber_nonfpt, " << num_jumps_fiber_nonfpt << endl; stats_results_per_file_number_fout << "total_particle_end_of_run, " << total_particle_end_of_run << ", "; stats_results_per_file_number_fout << "max_distance_exceeded_count, " << max_distance_exceeded_count; stats_results_per_file_number_fout << ",%, " << (float)max_distance_exceeded_count/(float)total_particle_end_of_run << endl; stats_results_per_file_number_fout << ",,displ_diag_exceeded_count, " << total_particle_end_of_run-max_distance_exceeded_count; stats_results_per_file_number_fout << ",%, " << (float)(total_particle_end_of_run-max_distance_exceeded_count)/(float)total_particle_end_of_run << endl; stats_results_per_file_number_fout << endl << "***results - buckets" << endl; stats_results_per_file_number_fout << "dist, dispx, dispy, dispz, %" << endl; // - least square fit - linear regression //http://mathworld.wolfram.com/LeastSquaresFitting.html //STEP1 //get bucket width: we want the buckets to cover the distance //betweeen 0 and a limited_max_distance, such that percent_max_distance of the datapoints are < limited_max_distance (and > 0 of course) //get limited_max_distance first: use a total_num_buckets buckets histogram to determine it // Grab start of list pointer LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); //get local values from dialog unsigned int total_num_buckets, bucket_index; total_num_buckets = dlg->field_rw_NUM_BUCKETS1_int; unsigned long long total_count_distance_bucket = 0, total_count_distance_bucket2 = 0;; float percent_max_distance; percent_max_distance = dlg->field_rw_PERCENT_MAX_DISTANCE_float; float percent_densest_bucket; percent_densest_bucket = dlg->field_rw_PERCENT_DENSEST_BUCKET_float; float limited_max_distance; float bucket_width; bucket_width = (max_distance+(float)1)/(float)total_num_buckets; //in microns unsigned long long * count_distance_bucket1 = new unsigned long long [total_num_buckets+1]; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket1[bucket_index] = (unsigned long long)0; } while ( distance_structure_list_node_ptr ) { distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket1[bucket_index] ++; //find max displ2's //max_displ2x = max(max_displ2x, temp_entry.dispx); //max_displ2y = max(max_displ2y, temp_entry.dispy); //max_displ2z = max(max_displ2z, temp_entry.dispz); distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } total_count_distance_bucket = 0; stats_results_per_file_number_fout << "raw data histogram, " << endl; stats_results_per_file_number_fout << "distance, bucket_count, " << endl; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { stats_results_per_file_number_fout << ((float)bucket_index+.5)*bucket_width << ", " << count_distance_bucket1[bucket_index] << ", " << endl; total_count_distance_bucket += count_distance_bucket1[bucket_index]; } if (dlg->stop_request) { if (total_count_distance_bucket < 200000) //500000 goto dostats;//delvar; //too few data points, quit now else rod_zero_aborted_enough_data_for_stats = true; //enough data, can proceed to doing stats } total_count_distance_bucket2 = 0; bucket_index=0; while (total_count_distance_bucket2 < percent_max_distance*total_count_distance_bucket) {//cumulative distribution of datapoints with increasing distance until criterion is reached ("percent_max_distance") total_count_distance_bucket2 += count_distance_bucket1[bucket_index]; bucket_index ++; } delete [] count_distance_bucket1; stats_results_per_file_number_fout << "total_bucket_count, " << total_count_distance_bucket << endl; stats_results_per_file_number_fout << "max_distance, " << max_distance << endl; limited_max_distance = ((float)bucket_index+1)*bucket_width; stats_results_per_file_number_fout << 100*percent_max_distance << "%_max_distance, " << limited_max_distance << endl; stats_results_per_file_number_fout << "total_particle_end_of_run, " << total_particle_end_of_run << ", "; stats_results_per_file_number_fout << "max_distance_exceeded_count, " << max_distance_exceeded_count; stats_results_per_file_number_fout << ",%, " << (float)max_distance_exceeded_count/(float)total_particle_end_of_run << endl; stats_results_per_file_number_fout << ",,displ_diag_exceeded_count, " << total_particle_end_of_run-max_distance_exceeded_count; stats_results_per_file_number_fout << ",%, " << (float)(total_particle_end_of_run-max_distance_exceeded_count)/(float)total_particle_end_of_run << endl; //STEP2 total_num_buckets = dlg->field_rw_NUM_BUCKETS2_int; unsigned long long * count_distance_bucket = new unsigned long long [total_num_buckets+1]; double * distance_bucket = new double [total_num_buckets+1]; double * disp2x_bucket = new double [total_num_buckets+1]; double * disp2y_bucket = new double [total_num_buckets+1]; double * disp2z_bucket = new double [total_num_buckets+1]; double * sd_distance_bucket = new double [total_num_buckets+1]; double * sd_dispx_bucket = new double [total_num_buckets+1]; double * sd_dispy_bucket = new double [total_num_buckets+1]; double * sd_dispz_bucket = new double [total_num_buckets+1]; unsigned short int num_used_buckets; double SSxx_distance, SSyy_dispx, SSyy_dispy, SSyy_dispz; //(SSxx and SSyys) unsigned long long count_distance_bucket_for_grand_av; double SSxx_distance_for_grand_av; double SSxy_ddx, SSxy_ddy, SSxy_ddz; //(SSxys) float av_distance, av_dispx, av_dispy, av_dispz; do { stats_results_per_file_number_fout << "***results - buckets, " << 100*percent_densest_bucket << ", % of densest bucket, " << endl; bucket_width = limited_max_distance/(float)total_num_buckets; float num_times_dataCollectionPoint_increment; num_times_dataCollectionPoint_increment = (float)1; if (bucket_width < num_times_dataCollectionPoint_increment*dataCollectionPoint_increment) { write_to_log("bucket_width @ STEP2 too small = %f", bucket_width); bucket_width = num_times_dataCollectionPoint_increment*dataCollectionPoint_increment; //to ensure enough write_to_log("increasing size to %f *dataCollectionPoint_increment = %f", num_times_dataCollectionPoint_increment, bucket_width); } stats_results_per_file_number_fout << "max_num_dist_buckets, " << total_num_buckets << ", dist_bucket_width, " << bucket_width << endl; //write these lines later //stats_results_per_file_number_fout << "(dist in microns), (displ^2 in microns^2)" << endl; //stats_results_per_file_number_fout << "av_dist, av_displ2X, av_displ2Y, av_displ2Z, count, "; //stats_results_per_file_number_fout << "sd_dist_norm'd, sd_disp2X_norm'd, sd_disp2Y_norm'd, sd_disp2Z_norm'd, "; ////stats_results_per_file_number_fout << "(av-sd)_dist, (av+sd)_dist, (av-sd)_disp2X, (av+sd)_disp2X, "; ////stats_results_per_file_number_fout << "(av-sd)_disp2Y, (av+sd)_disp2Y, (av-sd)_disp2Z, (av+sd)_disp2Z, "; //stats_results_per_file_number_fout << "k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; stats_results_per_file_number_fout << endl; unsigned long long max_num_in_a_bucket = 0; //unsigned long long counter; //count_distances //total_num_buckets = (unsigned short int)(limited_max_distance/bucket_width) + 1; //float av_buckets; //unsigned short int num_av_buckets; //av_buckets = (float)0; //num_av_buckets = 0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket[bucket_index] = (unsigned long long)0; distance_bucket[bucket_index] = (double)0; disp2x_bucket[bucket_index] = (double)0; disp2y_bucket[bucket_index] = (double)0; disp2z_bucket[bucket_index] = (double)0; sd_distance_bucket[bucket_index] = (double)0; sd_dispx_bucket[bucket_index] = (double)0; sd_dispy_bucket[bucket_index] = (double)0; sd_dispz_bucket[bucket_index] = (double)0; } //LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); if (save_full_RW_data) //will save the RAW data, unfiltered (if bucket counts are too small they will still be present here) { //create array to store displ2_xyz in buckets of fixed distance (simple histogram) DS_ull_cube * ds_cube_dist_displ2; int total_num_buckets_displ2;// = 500; //500 * 100 = 50 000 < 65535 -> csv file can be opened in excel (ms office <= 2003) total_num_buckets_displ2 = (int)((double)50000 / (double)total_num_buckets); //more displ^2 buckets if fewer dist buckets for greater detail //for simplicity - can be improved for better precision //max_displ2 = max_displacement_diag_pixels_sq; //should probably be chosen smaller for better precision max_displ2 = sq_MAX_ROI_RADIUS * PIXEL_SIZE_sq / 2; //in microns^2 // /2 to get better accuracy float buckets_displ2_width = max_displ2 / (float)total_num_buckets_displ2; float buckets_displ2_width_inv = (float)total_num_buckets_displ2 / max_displ2; //max_displ2x = max_displ2y = max_displ2z = max( max_displ2x, max(max_displ2y, max_displ2z) ); //float buckets_displ2x_width = max_displ2x / (float)total_num_buckets_displ2; //float buckets_displ2x_width_inv = (float)total_num_buckets_displ2 / max_displ2x; //float buckets_displ2y_width = max_displ2y / (float)total_num_buckets_displ2; //float buckets_displ2y_width_inv = (float)total_num_buckets_displ2 / max_displ2y; //float buckets_displ2z_width = max_displ2z / (float)total_num_buckets_displ2; //float buckets_displ2z_width_inv = (float)total_num_buckets_displ2 / max_displ2z; ds_cube_dist_displ2 = new DS_ull_cube( 3, (unsigned short)(total_num_buckets+1), (unsigned short)(total_num_buckets_displ2+1) ); //warning: no more than 65000 displ2 buckets! (unsigned short) //4=dist, 3=displ2x,displ2y,displ2z //total_num_buckets # of dist buckets //total_num_buckets_displ2 # of displ2 buckets ds_cube_dist_displ2->clear(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc bucket AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); //dist bucket count_distance_bucket[bucket_index] ++; distance_bucket[bucket_index] += (double)temp_entry.distance; disp2x_bucket[bucket_index] += (double)temp_entry.dispx; disp2y_bucket[bucket_index] += (double)temp_entry.dispy; disp2z_bucket[bucket_index] += (double)temp_entry.dispz; //add ins to save the displ2 histogram ds_cube_dist_displ2->data[0][bucket_index][(int)(temp_entry.dispx*buckets_displ2_width_inv)] ++; ds_cube_dist_displ2->data[1][bucket_index][(int)(temp_entry.dispy*buckets_displ2_width_inv)] ++; ds_cube_dist_displ2->data[2][bucket_index][(int)(temp_entry.dispz*buckets_displ2_width_inv)] ++; //ds_cube_dist_displ2->data[0][bucket_index][(int)(temp_entry.dispx*buckets_displ2x_width_inv)] ++; //ds_cube_dist_displ2->data[1][bucket_index][(int)(temp_entry.dispy*buckets_displ2y_width_inv)] ++; //ds_cube_dist_displ2->data[2][bucket_index][(int)(temp_entry.dispz*buckets_displ2z_width_inv)] ++; //ds_cube_dist_displ2[3][bucket_index][] ++; } //no save_full_RW_data, next datapoint in the distance structure distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } char rw_raw_data_filename[1024]; char* rw_raw_data_file; //sprintf(rw_raw_data_filename, "%s_rw_raw_data_bucket#%d.csv", output_full_name, bucket_index); sprintf(rw_raw_data_filename, "%s_rw_raw_data_per_dist_bucket.csv", output_full_name); rw_raw_data_file = rw_raw_data_filename; rw_raw_data_fout.open(rw_raw_data_file, ios::out | ios::app); if( !rw_raw_data_fout.is_open() ) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving rw_raw_data: %s", rw_raw_data_filename); write_to_log(text); return; } //file header rw_raw_data_fout << output_full_name << endl; //sample name rw_raw_data_fout << "bucket # of displ2_, av_displ2_, _x#, _y#, _z#, gammaX, gammaY, gammaZ"; //count-histogram rw_raw_data_fout << endl; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { //file sub-header / per dist bucket rw_raw_data_fout << "dist_bucket#, " << bucket_index << ", av_dist, " << distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index] ; rw_raw_data_fout << endl; int i; double displ2_data0_av, displ2_data1_av, displ2_data2_av;//0=X, 1=Y, 2=Z double displ2_data0_var, displ2_data1_var, displ2_data2_var; unsigned long long count_num_datapoints = 0; { //calc displ2 bucket avg displ2_data0_av = displ2_data1_av = displ2_data2_av = 0; for (i=0; i<total_num_buckets_displ2; i++) { displ2_data0_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[0][bucket_index][i]; displ2_data1_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[1][bucket_index][i]; displ2_data2_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[2][bucket_index][i]; count_num_datapoints += ds_cube_dist_displ2->data[0][bucket_index][i]; //sum should be equal for x, y and z data } displ2_data0_av /= (double)count_num_datapoints; displ2_data1_av /= (double)count_num_datapoints; displ2_data2_av /= (double)count_num_datapoints; } //calc displ2 bucket variance displ2_data0_var = displ2_data1_var = displ2_data2_var = 0; for (i=0; i<total_num_buckets_displ2; i++) { displ2_data0_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[0][bucket_index][i] - displ2_data0_av, 2); displ2_data1_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[1][bucket_index][i] - displ2_data1_av, 2); displ2_data2_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[2][bucket_index][i] - displ2_data2_av, 2); } displ2_data0_av *= buckets_displ2_width; displ2_data1_av *= buckets_displ2_width; displ2_data2_av *= buckets_displ2_width; displ2_data0_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; displ2_data1_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; displ2_data2_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; rw_raw_data_fout << "displ2 av_, _x, _y, _z" << endl; rw_raw_data_fout << ", " << displ2_data0_av; rw_raw_data_fout << ", " << displ2_data1_av; rw_raw_data_fout << ", " << displ2_data2_av; rw_raw_data_fout << endl; rw_raw_data_fout << "displ2 sd_, _x, _y, _z" << endl; rw_raw_data_fout << ", " << sqrt(displ2_data0_var); rw_raw_data_fout << ", " << sqrt(displ2_data1_var); rw_raw_data_fout << ", " << sqrt(displ2_data2_var); rw_raw_data_fout << endl; //estimate gamma distrib parameters double alphax, betax, alphay, betay, alphaz, betaz; alphax = pow(displ2_data0_av, 2) / displ2_data0_var; betax = displ2_data0_var / displ2_data0_av; alphay = pow(displ2_data1_av, 2) / displ2_data1_var; betay = displ2_data1_var / displ2_data1_av; alphaz = pow(displ2_data2_av, 2) / displ2_data2_var; betaz = displ2_data2_var / displ2_data2_av; rw_raw_data_fout << ", k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; rw_raw_data_fout << endl; rw_raw_data_fout << ", " << alphax << ", " << betax; rw_raw_data_fout << ", " << alphay << ", " << betay; rw_raw_data_fout << ", " << alphaz << ", " << betaz; rw_raw_data_fout << endl; //save whole histogram data to file + add Gamma distribution values for (i=0; i<total_num_buckets_displ2; i++) { double t = (double)(i+.5)*(double)buckets_displ2_width; rw_raw_data_fout << i << ","; rw_raw_data_fout << t << ","; //middle of displ2 bucket av //rw_raw_data_fout << << ","; //true av for displ2_x is not available here - would require extra structure to hold data rw_raw_data_fout << ds_cube_dist_displ2->data[0][bucket_index][i] << ","; rw_raw_data_fout << ds_cube_dist_displ2->data[1][bucket_index][i] << ","; rw_raw_data_fout << ds_cube_dist_displ2->data[2][bucket_index][i] << ","; rw_raw_data_fout << gamma_pdf(t, alphax, betax) << ","; rw_raw_data_fout << gamma_pdf(t, alphay, betay) << ","; rw_raw_data_fout << gamma_pdf(t, alphaz, betaz) << ","; rw_raw_data_fout << endl; } //file sub-footer / per dist bucket rw_raw_data_fout << "********************************************************************************************************"; rw_raw_data_fout << endl << endl; } //close raw rw file rw_raw_data_fout << endl;// << "****" << endl << endl; rw_raw_data_fout.close(); //delete histogram data cube delete ds_cube_dist_displ2; } else { while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket[bucket_index] ++; distance_bucket[bucket_index] += (double)temp_entry.distance; disp2x_bucket[bucket_index] += (double)temp_entry.dispx; disp2y_bucket[bucket_index] += (double)temp_entry.dispy; disp2z_bucket[bucket_index] += (double)temp_entry.dispz; } //no save_full_RW_data, next datapoint in the distance structure distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } } //while ( distance_structure_list_node_ptr ) //{//summing everything - _bucket is just the SUM here -> to calc AVERAGES first // distance_structure temp_entry; // temp_entry = distance_structure_list_node_ptr->item; // if (temp_entry.distance < limited_max_distance) // { // bucket_index = (unsigned int)(temp_entry.distance/bucket_width); // count_distance_bucket[bucket_index] ++; // distance_bucket[bucket_index] += (double)temp_entry.distance; // disp2x_bucket[bucket_index] += (double)temp_entry.dispx; // disp2y_bucket[bucket_index] += (double)temp_entry.dispy; // disp2z_bucket[bucket_index] += (double)temp_entry.dispz; // } // distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; //} //// should depend on total_num_buckets? // depending on a fixed value -> not good //min_required_num_per_bucket = (unsigned int)(count_distances*0.02); //xx% of total # of datapoints -> not good either // depending on value in most populated bucket for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { if ( max_num_in_a_bucket < count_distance_bucket[bucket_index] ) max_num_in_a_bucket = count_distance_bucket[bucket_index]; } min_required_num_per_bucket = (unsigned int)(max_num_in_a_bucket*percent_densest_bucket); //percent_densest_bucket% of most populated bucket num_used_buckets = 0; count_distance_bucket_for_grand_av = 0; total_count_distance_bucket = 0; SSxx_distance = SSxx_distance_for_grand_av = (double)0; SSyy_dispx = SSyy_dispy = SSyy_dispz = (double)0; SSxy_ddx = SSxy_ddy = SSxy_ddz = (double)0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//summing everything - SS is just the overall SUM here -> temp variables to calc overall AVERAGES first if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) {//keep value only if there are enough datapoints in that particular bucket //this is to calculate the grand_av - expected to be much smaller than the lin reg av SSxx_distance_for_grand_av += distance_bucket[bucket_index]; SSxy_ddx += disp2x_bucket[bucket_index]; SSxy_ddy += disp2y_bucket[bucket_index]; //these are not SSxys they are temp variables reused to calc the grand_av SSxy_ddz += disp2z_bucket[bucket_index]; count_distance_bucket_for_grand_av += count_distance_bucket[bucket_index]; //SSxx_distance += distance_bucket[bucket_index] = ((double)bucket_index+.5) * (double)bucket_width; //SSxx_distance calculation shouldn't involve middle of bucket value, instead better to use real data value: SSxx_distance += distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //this is calculating avrg for lin regression of ensemble averages (ie per dist bucket - not grand_av) SSyy_dispx += disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispy += disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispz += disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //av per bucket - calc'd above already //distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; num_used_buckets ++; total_count_distance_bucket += count_distance_bucket[bucket_index]; } else { distance_bucket[bucket_index] = 0; disp2x_bucket[bucket_index] = 0; //averages per bucket disp2y_bucket[bucket_index] = 0; disp2z_bucket[bucket_index] = 0; } } percent_densest_bucket -= (float).01; //decrease % to ensure at least 3 buckets } while ( (num_used_buckets < 3) && (percent_densest_bucket >= 0) ); if (num_used_buckets == 1 ) { write_to_log("only 1 bucket!"); } //calc grand (overall) averages av_distance = (float)(SSxx_distance_for_grand_av/(double)count_distance_bucket_for_grand_av); //grand xbar av_dispx = (float)(SSxy_ddx/(double)count_distance_bucket_for_grand_av); av_dispy = (float)(SSxy_ddy/(double)count_distance_bucket_for_grand_av); //grand ybar av_dispz = (float)(SSxy_ddz/(double)count_distance_bucket_for_grand_av); stats_results_per_file_number_fout << "grand avg, " << endl; stats_results_per_file_number_fout << "dist, displ2X, displ2Y, displ2Z, " << endl; stats_results_per_file_number_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; //estimated D/D0 coefficients stats_results_per_file_number_fout << "and estimated D/D0 coefficients (xyz)" << endl; { D_D0x = (float)1.5*av_dispx/av_distance/lambda_micron; D_D0y = (float)1.5*av_dispy/av_distance/lambda_micron; D_D0z = (float)1.5*av_dispz/av_distance/lambda_micron; } stats_results_per_file_number_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << endl; stats_results_per_file_number_fout << endl; //calc lin reg avg av_distance = (float)SSxx_distance/(float)num_used_buckets; //xbar av_dispx = (float)SSyy_dispx/(float)num_used_buckets; av_dispy = (float)SSyy_dispy/(float)num_used_buckets; //ybar av_dispz = (float)SSyy_dispz/(float)num_used_buckets; stats_results_per_file_number_fout << "lin reg avg, " << endl; stats_results_per_file_number_fout << "dist, displ2X, displ2Y, displ2Z, " << endl; stats_results_per_file_number_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; //estimated D/D0 coefficients stats_results_per_file_number_fout << "and estimated D/D0 coefficients (xyz)" << endl; { D_D0x = (float)1.5*av_dispx/av_distance/lambda_micron; D_D0y = (float)1.5*av_dispy/av_distance/lambda_micron; D_D0z = (float)1.5*av_dispz/av_distance/lambda_micron; } stats_results_per_file_number_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << endl; stats_results_per_file_number_fout << endl; //now that we have the averages //calc standard errors for everything // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) { sd_distance_bucket[bucket_index] += pow((double)temp_entry.distance - distance_bucket[bucket_index], 2); //this is not the sd yet (see final step below) sd_dispx_bucket[bucket_index] += pow((double)temp_entry.dispx - disp2x_bucket[bucket_index], 2); sd_dispy_bucket[bucket_index] += pow((double)temp_entry.dispy - disp2y_bucket[bucket_index], 2); sd_dispz_bucket[bucket_index] += pow((double)temp_entry.dispz - disp2z_bucket[bucket_index], 2); } } distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } //sub part header stats_results_per_file_number_fout << "(dist in microns), (displ^2 in microns^2)" << endl; stats_results_per_file_number_fout << "av_dist, av_displ2X, av_displ2Y, av_displ2Z, count, "; stats_results_per_file_number_fout << "sd_dist_norm'd, sd_disp2X_norm'd, sd_disp2Y_norm'd, sd_disp2Z_norm'd, "; //stats_results_per_file_number_fout << "(av-sd)_dist, (av+sd)_dist, (av-sd)_disp2X, (av+sd)_disp2X, "; //stats_results_per_file_number_fout << "(av-sd)_disp2Y, (av+sd)_disp2Y, (av-sd)_disp2Z, (av+sd)_disp2Z, "; //stats_results_per_file_number_fout << "k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; stats_results_per_file_number_fout << endl; //init variables SSxx_distance = (double)0; //SSxx SSyy_dispx = (double)0; SSyy_dispy = (double)0; //SSyy SSyy_dispz = (double)0; SSxy_ddx = (double)0; SSxy_ddy = (double)0; //SSxy SSxy_ddz = (double)0; //calc SS for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) if (distance_bucket[bucket_index] > 0) { SSxx_distance += pow(distance_bucket[bucket_index], 2); //SSxx=SUM(x^2) SSyy_dispx += pow(disp2x_bucket[bucket_index], 2); SSyy_dispy += pow(disp2y_bucket[bucket_index], 2); //SSyy=SUM(y^2) SSyy_dispz += pow(disp2z_bucket[bucket_index], 2); SSxy_ddx += distance_bucket[bucket_index] * disp2x_bucket[bucket_index]; SSxy_ddy += distance_bucket[bucket_index] * disp2y_bucket[bucket_index]; //SSxy=SUM(x*y) SSxy_ddz += distance_bucket[bucket_index] * disp2z_bucket[bucket_index]; sd_distance_bucket[bucket_index] = sqrt(sd_distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispx_bucket[bucket_index] = sqrt(sd_dispx_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); //sd = sqrt(sum of var)/n sd_dispy_bucket[bucket_index] = sqrt(sd_dispy_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispz_bucket[bucket_index] = sqrt(sd_dispz_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); //write values to file - to check the lin. regression fit //"av_dist, av_dispx, av_dispy, av_dispz, count, sd_dist_norm, sd_dispx_norm, sd_dispy_norm, sd_dispz_norm, //(av-sd)_dist, (av+sd)_dist, (av-sd)_dispX, (av+sd)_dispX, (av-sd)_dispY, (av+sd)_dispY, (av-sd)_dispZ, (av+sd)_dispZ, k=alpha, theta=beta," stats_results_per_file_number_fout \ << distance_bucket[bucket_index] << ", " \ << disp2x_bucket[bucket_index] << ", " \ << disp2y_bucket[bucket_index] << ", " \ << disp2z_bucket[bucket_index] << ", " \ << count_distance_bucket[bucket_index] << ", "; stats_results_per_file_number_fout \ << (sd_distance_bucket[bucket_index]/distance_bucket[bucket_index]) << ", " \ << (sd_dispx_bucket[bucket_index]/disp2x_bucket[bucket_index]) << ", " \ << (sd_dispy_bucket[bucket_index]/disp2y_bucket[bucket_index]) << ", " \ << (sd_dispz_bucket[bucket_index]/disp2z_bucket[bucket_index]) << ", "; //stats_results_per_file_number_fout \ // << (distance_bucket[bucket_index] - sd_distance_bucket[bucket_index]) << ", " \ // << (distance_bucket[bucket_index] + sd_distance_bucket[bucket_index]) << ", " \ // << (disp2x_bucket[bucket_index] - sd_dispx_bucket[bucket_index]) << ", " \ // << (disp2x_bucket[bucket_index] + sd_dispx_bucket[bucket_index]) << ", " \ // << (disp2y_bucket[bucket_index] - sd_dispy_bucket[bucket_index]) << ", " \ // << (disp2y_bucket[bucket_index] + sd_dispy_bucket[bucket_index]) << ", " \ // << (disp2z_bucket[bucket_index] - sd_dispz_bucket[bucket_index]) << ", " \ // << (disp2z_bucket[bucket_index] + sd_dispz_bucket[bucket_index]) << ", "; stats_results_per_file_number_fout << endl; } delete [] count_distance_bucket; delete [] distance_bucket; delete [] disp2x_bucket; delete [] disp2y_bucket; delete [] disp2z_bucket; delete [] sd_distance_bucket; delete [] sd_dispx_bucket; delete [] sd_dispy_bucket; delete [] sd_dispz_bucket; //save results for 0-intercept float Sumxx, Sumxyx, Sumxyy, Sumxyz; Sumxx = (float)SSxx_distance; Sumxyx = (float)SSxy_ddx; Sumxyy = (float)SSxy_ddy; Sumxyz = (float)SSxy_ddz; //finally: SSxx_distance = SSxx_distance - (double)num_used_buckets * pow((double)av_distance, 2); //SSxx=SUM(x^2)-n*xbar^2 SSyy_dispx = SSyy_dispx - (double)num_used_buckets * pow((double)av_dispx, 2); SSyy_dispy = SSyy_dispy - (double)num_used_buckets * pow((double)av_dispy, 2); //SSyy=SUM(y^2)-n*ybar^2 SSyy_dispz = SSyy_dispz - (double)num_used_buckets * pow((double)av_dispz, 2); SSxy_ddx = SSxy_ddx - (double)num_used_buckets * (double)av_distance * (double)av_dispx; SSxy_ddy = SSxy_ddy - (double)num_used_buckets * (double)av_distance * (double)av_dispy; //SSxy=SUM(x*y)-n*xbar*ybar SSxy_ddz = SSxy_ddz - (double)num_used_buckets * (double)av_distance * (double)av_dispz; //get a, b (y=a+bx) -- disp^2 vs distance //b=SSxy/SSxx float bx, by, bz; bx = (float)(SSxy_ddx/SSxx_distance); by = (float)(SSxy_ddy/SSxx_distance); bz = (float)(SSxy_ddz/SSxx_distance); //(all we want is the slope b) //now a=ybar-b*xbar float ax, ay, az; ax = av_dispx - bx*av_distance; ay = av_dispy - by*av_distance; az = av_dispz - bz*av_distance; float sigmax_distance, sigmay_dispx, sigmay_dispy, sigmay_dispz, covxy_ddx, covxy_ddy, covxy_ddz; //sigmax^2=SSxx/n sigmax_distance = (float)(SSxx_distance / num_used_buckets); //sigmay^2=SSyy/n sigmay_dispx = (float)(SSyy_dispx / (double)num_used_buckets); sigmay_dispy = (float)(SSyy_dispy / (double)num_used_buckets); sigmay_dispz = (float)(SSyy_dispz / (double)num_used_buckets); //Cov(x,y)=SSxy/n covxy_ddx = (float)(SSxy_ddx / (double)num_used_buckets); covxy_ddy = (float)(SSxy_ddy / (double)num_used_buckets); covxy_ddz = (float)(SSxy_ddz / (double)num_used_buckets); //s=sqrt( (SSyy-b*SSxy)/(n-2) ) float s_x, s_y, s_z; s_x = (float)sqrt( (SSyy_dispx-bx*SSxy_ddx)/((float)num_used_buckets-2) ); s_y = (float)sqrt( (SSyy_dispy-by*SSxy_ddy)/((float)num_used_buckets-2) ); s_z = (float)sqrt( (SSyy_dispz-bz*SSxy_ddz)/((float)num_used_buckets-2) ); //se(a)=s*sqrt( 1/n + xbar^2/SSxx ) float se_ax, se_ay, se_az; se_ax = s_x*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_ay = s_y*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_az = s_z*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); //se(b)=s/sqrt(SSxx) float se_bx, se_by, se_bz; se_bx = s_x / (float)sqrt(SSxx_distance); se_by = s_y / (float)sqrt(SSxx_distance); se_bz = s_z / (float)sqrt(SSxx_distance); //R^2=SSxy^2/SSxx/SSyy float R2_x, R2_y, R2_z; R2_x = (float)(SSxy_ddx*SSxy_ddx/SSxx_distance/SSyy_dispx); R2_y = (float)(SSxy_ddy*SSxy_ddy/SSxx_distance/SSyy_dispy); R2_z = (float)(SSxy_ddz*SSxy_ddz/SSxx_distance/SSyy_dispz); //store stats results to file stats_results_per_file_number_fout << "average" << endl; stats_results_per_file_number_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; stats_results_per_file_number_fout << "sigma^2" << endl; stats_results_per_file_number_fout << sigmax_distance << ", " << sigmay_dispx << ", " << sigmay_dispy << ", " << sigmay_dispz << endl << endl; stats_results_per_file_number_fout << "(y=a+bx) -- disp^2 vs distance" << endl; stats_results_per_file_number_fout << ", disp2x, disp2y, disp2z, " << endl; stats_results_per_file_number_fout << "a, " << ax << ", " << ay << ", " << az << endl; stats_results_per_file_number_fout << "b, " << bx << ", " << by << ", " << bz << endl; stats_results_per_file_number_fout << "R^2, " << R2_x << ", " << R2_y << ", " << R2_z << endl; stats_results_per_file_number_fout << "se_a, " << se_ax << ", " << se_ay << ", " << se_az << endl; stats_results_per_file_number_fout << "se_b, " << se_bx << ", " << se_by << ", " << se_bz << endl; //null intercept, b=Sumxy[]/Sumxx stats_results_per_file_number_fout << "(y=bx) -- zero (0) intercept" << endl; stats_results_per_file_number_fout << "b, " << (float)(Sumxyx/Sumxx) \ << ", " << (float)(Sumxyy/Sumxx) \ << ", " << (float)(Sumxyz/Sumxx) << endl; //now get D/D0 coefs stats_results_per_file_number_fout << endl; //stats_results_per_file_number_fout << "D/D0 coefficients (xyz),,,,,coeff/porosity" << endl; stats_results_per_file_number_fout << "D/D0 coefficients (xyz)," << endl; { D_D0x = (float)1.5*bx*ROI_porosity/lambda_micron; D_D0y = (float)1.5*by*ROI_porosity/lambda_micron; D_D0z = (float)1.5*bz*ROI_porosity/lambda_micron; } stats_results_per_file_number_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << endl; //stats_results_per_file_number_fout << ", " << D_D0x/ROI_porosity << ", " << D_D0y/ROI_porosity << ", " << D_D0z/ROI_porosity << endl; stats_results_per_file_number_fout << "+/-" << endl; stats_results_per_file_number_fout \ << ", " << (float)1.5*se_bx*ROI_porosity/lambda_micron \ << ", " << (float)1.5*se_by*ROI_porosity/lambda_micron \ << ", " << (float)1.5*se_bz*ROI_porosity/lambda_micron << endl; stats_results_per_file_number_fout << "with zero (0) intercept" << endl; stats_results_per_file_number_fout \ << ", " << (float)1.5*ROI_porosity*(Sumxyx/Sumxx)/lambda_micron \ << ", " << (float)1.5*ROI_porosity*(Sumxyy/Sumxx)/lambda_micron \ << ", " << (float)1.5*ROI_porosity*(Sumxyz/Sumxx)/lambda_micron << endl; stats_results_per_file_number_fout << endl; {//sum up coefs for final result (average in stats_results_fout) coef_rod_x += (float)1.5*bx*ROI_porosity/lambda_micron; coef_rod_y += (float)1.5*by*ROI_porosity/lambda_micron; coef_rod_z += (float)1.5*bz*ROI_porosity/lambda_micron; num_coeff ++; coef_rod_sex += (float)1.5*se_bx*ROI_porosity/lambda_micron; coef_rod_sey += (float)1.5*se_by*ROI_porosity/lambda_micron; coef_rod_sez += (float)1.5*se_bz*ROI_porosity/lambda_micron; coef_rod0_x = (float)1.5*ROI_porosity*(Sumxyx/Sumxx)/lambda_micron; coef_rod0_y = (float)1.5*ROI_porosity*(Sumxyy/Sumxx)/lambda_micron; coef_rod0_z = (float)1.5*ROI_porosity*(Sumxyz/Sumxx)/lambda_micron; } stats_results_per_file_number_fout << endl << "****" << endl << endl; stats_results_per_file_number_fout.close(); //clear all variables that need to be reused distance_structure_list_ptr->removeAllNodes(); // distance_structure_list_ptr->purgeFreeList(); //because we dont need them anymore afterwards max_distance = (float)0; count_distances = 0; } if (dlg->stop_request) { rw_cv_data_fout << "CANCELED @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_fout.close(); //goto delvar; total_count_tracer_start_hit_pore += count_tracer_start_hit_pore; total_count_tracer_start_hit += count_tracer_start_hit; goto dostats; } total_count_tracer_start_hit_pore += count_tracer_start_hit_pore; total_count_tracer_start_hit += count_tracer_start_hit; } dostats: if (dlg->stop_request) { sprintf(win_txt, "stop_request received, attempting to do stats and finish now"); dlg->txt_field_progress.SetWindowText(win_txt); write_to_log(win_txt); } write_to_log("total_particle_end_of_run, %d", total_particle_end_of_run); //do stats if (rod) //ie != 0 {//do stats if (!dlg->stop_request) { write_to_log("Simulation: done - now doing final stats"); dlg->txt_field_progress.SetWindowText("Simulation: done - now doing stats"); } float average_ROI_porosity = total_ROI_porosity/(float)total_file_number; stats_results_fout << "average_true_ROI_porosity, " << average_ROI_porosity << endl; stats_results_fout << "average_ROI_porosity_from tracer start, " \ << (float)total_count_tracer_start_hit_pore/(float)total_count_tracer_start_hit << endl; // - least square fit - linear regression //http://mathworld.wolfram.com/LeastSquaresFitting.html //STEP1 //get bucket width: we want the buckets to cover the distance //betweeen 0 and a limited_max_distance, such that percent_max_distance of the datapoints are < limited_max_distance (and > 0 of course) //get limited_max_distance first: use a total_num_buckets buckets histogram to determine it // Grab start of list pointer LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); unsigned int total_num_buckets, bucket_index; total_num_buckets = dlg->field_rw_NUM_BUCKETS1_int;//10000; unsigned long long total_count_distance_bucket = 0, total_count_distance_bucket2 = 0; float percent_max_distance; percent_max_distance = dlg->field_rw_PERCENT_MAX_DISTANCE_float;//(float).5; float percent_densest_bucket; percent_densest_bucket = dlg->field_rw_PERCENT_DENSEST_BUCKET_float;//(float).9; float limited_max_distance; float bucket_width; bucket_width = (max_distance+(float)1)/(float)total_num_buckets; //in microns unsigned long long * count_distance_bucket1 = new unsigned long long [total_num_buckets+1]; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket1[bucket_index] = (unsigned long long)0; } while ( distance_structure_list_node_ptr ) { distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket1[bucket_index] ++; distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } total_count_distance_bucket = 0; stats_results_fout << "raw data histogram, " << endl; stats_results_fout << "distance, bucket_count, " << endl; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { stats_results_fout << ((float)bucket_index+.5)*bucket_width << ", " << count_distance_bucket1[bucket_index] << ", " << endl; total_count_distance_bucket += count_distance_bucket1[bucket_index]; } if (dlg->stop_request) { if (total_count_distance_bucket < 500000) goto delvar; //too few data points, quit now } total_count_distance_bucket2 = 0; bucket_index=0; while (total_count_distance_bucket2 < percent_max_distance*total_count_distance_bucket) { total_count_distance_bucket2 += count_distance_bucket1[bucket_index]; bucket_index ++; } delete [] count_distance_bucket1; stats_results_fout << "total_bucket_count, " << total_count_distance_bucket << endl; stats_results_fout << "max_distance, " << max_distance << endl; limited_max_distance = ((float)bucket_index+1)*bucket_width; stats_results_fout << 100*percent_max_distance << "%_max_distance, " << limited_max_distance << endl; stats_results_fout << "total_particle_end_of_run, " << total_particle_end_of_run << ", "; stats_results_fout << "max_distance_exceeded_count, " << max_distance_exceeded_count; stats_results_fout << ",%, " << (float)max_distance_exceeded_count/(float)total_particle_end_of_run << endl; stats_results_fout << ",,displ_diag_exceeded_count, " << total_particle_end_of_run-max_distance_exceeded_count; stats_results_fout << ",%, " << (float)(total_particle_end_of_run-max_distance_exceeded_count)/(float)total_particle_end_of_run << endl; //STEP2 total_num_buckets = dlg->field_rw_NUM_BUCKETS2_int;//20; unsigned long long * count_distance_bucket = new unsigned long long [total_num_buckets+1]; double * distance_bucket = new double [total_num_buckets+1]; double * disp2x_bucket = new double [total_num_buckets+1]; double * disp2y_bucket = new double [total_num_buckets+1]; double * disp2z_bucket = new double [total_num_buckets+1]; double * sd_distance_bucket = new double [total_num_buckets+1]; double * sd_dispx_bucket = new double [total_num_buckets+1]; double * sd_dispy_bucket = new double [total_num_buckets+1]; double * sd_dispz_bucket = new double [total_num_buckets+1]; unsigned short int num_used_buckets; unsigned long long count_distance_bucket_for_grand_av; double SSxx_distance_for_grand_av; double SSxx_distance, SSyy_dispx, SSyy_dispy, SSyy_dispz; //(SSxx and SSyys) double SSxy_ddx, SSxy_ddy, SSxy_ddz; //(SSxys) float av_distance, av_dispx, av_dispy, av_dispz; do { stats_results_fout << "***results - buckets, " << 100*percent_densest_bucket << ", % of densest bucket, " << endl; bucket_width = limited_max_distance/(float)total_num_buckets; float num_times_dataCollectionPoint_increment; num_times_dataCollectionPoint_increment = (float)1; if (bucket_width < num_times_dataCollectionPoint_increment*dataCollectionPoint_increment) { write_to_log("bucket_width @ STEP2 too small = %f", bucket_width); bucket_width = num_times_dataCollectionPoint_increment*dataCollectionPoint_increment; //to ensure enough write_to_log("increasing size to %f *dataCollectionPoint_increment = %f", num_times_dataCollectionPoint_increment, bucket_width); } stats_results_fout << "max_num_buckets, " << total_num_buckets << ", bucket_width, " << bucket_width << endl; //stats_results_fout << "(dist in microns), (displ^2 in microns^2)" << endl; //stats_results_fout << "av_dist, av_displ2X, av_displ2Y, av_displ2Z, count, "; //stats_results_fout << "sd_dist_norm, sd_disp2X_norm, sd_disp2Y_norm, sd_disp2Z_norm, "; //stats_results_fout << "(av-sd)_dist, (av+sd)_dist, (av-sd)_disp2X, (av+sd)_disp2X, "; //non normal data, so useless and removed //stats_results_fout << "(av-sd)_disp2Y, (av+sd)_disp2Y, (av-sd)_disp2Z, (av+sd)_disp2Z,"; //stats_results_fout << "k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; //trying to fit a gamma distribution stats_results_fout << endl; unsigned long long max_num_in_a_bucket = 0; //unsigned long long counter; //count_distances //total_num_buckets = (unsigned short int)(limited_max_distance/bucket_width) + 1; //float av_buckets; //unsigned short int num_av_buckets; //av_buckets = (float)0; //num_av_buckets = 0; SSxx_distance = (double)0; SSyy_dispx = (double)0; SSyy_dispy = (double)0; SSyy_dispz = (double)0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise couting variables count_distance_bucket[bucket_index] = (unsigned long long)0; distance_bucket[bucket_index] = (double)0; disp2x_bucket[bucket_index] = (double)0; disp2y_bucket[bucket_index] = (double)0; disp2z_bucket[bucket_index] = (double)0; sd_distance_bucket[bucket_index] = (double)0; sd_dispx_bucket[bucket_index] = (double)0; sd_dispy_bucket[bucket_index] = (double)0; sd_dispz_bucket[bucket_index] = (double)0; } //LinkedListNode<distance_structure> *distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); //if (save_full_RW_data) //{ // while ( distance_structure_list_node_ptr ) // {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first // distance_structure temp_entry; // temp_entry = distance_structure_list_node_ptr->item; // if (temp_entry.distance < limited_max_distance) // { // bucket_index = (unsigned int)(temp_entry.distance/bucket_width); //round distance for bucket # // count_distance_bucket[bucket_index] ++; // distance_bucket[bucket_index] += (double)temp_entry.distance; // disp2x_bucket[bucket_index] += (double)temp_entry.dispx; // disp2y_bucket[bucket_index] += (double)temp_entry.dispy; // disp2z_bucket[bucket_index] += (double)temp_entry.dispz; // } // //save_full_RW_data // rw_raw_data_fout << temp_entry.distance << ","; // rw_raw_data_fout << temp_entry.dispx << ","; // rw_raw_data_fout << temp_entry.dispy << ","; // rw_raw_data_fout << temp_entry.dispz << ","; // rw_raw_data_fout << endl; // //next datapoint in the distance structure // distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; // } // //close raw rw file // rw_raw_data_fout << endl;// << "****" << endl << endl; // rw_raw_data_fout.close(); //} if (save_full_RW_data) //will save the RAW data, unfiltered (if bucket counts are too small they will still be present here) { //create array to store displ2_xyz in buckets of fixed distance (simple histogram) DS_ull_cube * ds_cube_dist_displ2; DS_f_cube * ds_cube_min_displ2_per_dist_bucket; int total_num_buckets_displ2;// = 1000; //< 65535 -> csv file can be opened in excel (ms office <= 2003) total_num_buckets_displ2 = (int)((double)50000 / (double)total_num_buckets); //more displ^2 buckets if fewer dist buckets for greater detail max_displ2 = sq_MAX_ROI_RADIUS * PIXEL_SIZE_sq / 2; //in microns^2 // /2 to get better accuracy float buckets_displ2_width = max_displ2 / (float)total_num_buckets_displ2; float buckets_displ2_width_inv = (float)total_num_buckets_displ2 / max_displ2; ds_cube_dist_displ2 = new DS_ull_cube( 3, (unsigned short)(total_num_buckets+1), (unsigned short)(total_num_buckets_displ2+1) ); //warning: no more than 65000 displ2 buckets! (unsigned short) ds_cube_min_displ2_per_dist_bucket = new DS_f_cube( 1, 3, (unsigned short)(total_num_buckets+1)); //3=displ2x,displ2y,displ2z //total_num_buckets # of dist buckets //total_num_buckets_displ2 # of displ2 buckets for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//initialise structure to save min displ2 in each distance bucket ds_cube_min_displ2_per_dist_bucket->data[0][0][bucket_index] = pow((float)10, (float)38); ds_cube_min_displ2_per_dist_bucket->data[0][1][bucket_index] = pow((float)10, (float)38); ds_cube_min_displ2_per_dist_bucket->data[0][2][bucket_index] = pow((float)10, (float)38); } ds_cube_dist_displ2->clear(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); //dist bucket count_distance_bucket[bucket_index] ++; distance_bucket[bucket_index] += (double)temp_entry.distance; disp2x_bucket[bucket_index] += (double)temp_entry.dispx; disp2y_bucket[bucket_index] += (double)temp_entry.dispy; disp2z_bucket[bucket_index] += (double)temp_entry.dispz; //save min displ2 value if (temp_entry.dispx < ds_cube_min_displ2_per_dist_bucket->data[0][0][bucket_index]) ds_cube_min_displ2_per_dist_bucket->data[0][0][bucket_index] = temp_entry.dispx; if (temp_entry.dispy < ds_cube_min_displ2_per_dist_bucket->data[0][1][bucket_index]) ds_cube_min_displ2_per_dist_bucket->data[0][1][bucket_index] = temp_entry.dispy; if (temp_entry.dispz < ds_cube_min_displ2_per_dist_bucket->data[0][2][bucket_index]) ds_cube_min_displ2_per_dist_bucket->data[0][2][bucket_index] = temp_entry.dispz; ds_cube_dist_displ2->data[0][bucket_index][(int)(temp_entry.dispx*buckets_displ2_width_inv)] ++; ds_cube_dist_displ2->data[1][bucket_index][(int)(temp_entry.dispy*buckets_displ2_width_inv)] ++; ds_cube_dist_displ2->data[2][bucket_index][(int)(temp_entry.dispz*buckets_displ2_width_inv)] ++; } distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } char rw_raw_data_filename[1024]; char* rw_raw_data_file; //sprintf(rw_raw_data_filename, "%s_rw_raw_data_bucket#%d.csv", output_full_name, bucket_index); sprintf(rw_raw_data_filename, "%s_rw_raw_data_per_dist_bucket.csv", output_full_name); rw_raw_data_file = rw_raw_data_filename; rw_raw_data_fout.open(rw_raw_data_file, ios::out | ios::app); if( !rw_raw_data_fout.is_open() ) { char text[1024]; sprintf(text,"Simultaneous Diffusion: error opening file for saving rw_raw_data: %s", rw_raw_data_filename); write_to_log(text); return; } //file header rw_raw_data_fout << output_full_name << endl; //sample name rw_raw_data_fout << "bucket # of displ2_, av_displ2_, _x#, _y#, _z#, gammaX, gammaY, gammaZ"; //count-histogram rw_raw_data_fout << endl; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) { //file sub-header / per dist bucket rw_raw_data_fout << "dist bucket #, " << bucket_index << ", av_dist, " << distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index] ; rw_raw_data_fout << endl; rw_raw_data_fout << ", _x, _y, _z" << endl; rw_raw_data_fout << "min_displ2,"; rw_raw_data_fout << ds_cube_min_displ2_per_dist_bucket->data[0][0][bucket_index] << ","; rw_raw_data_fout << ds_cube_min_displ2_per_dist_bucket->data[0][1][bucket_index] << ","; rw_raw_data_fout << ds_cube_min_displ2_per_dist_bucket->data[0][2][bucket_index] << ","; rw_raw_data_fout << endl; int i; double displ2_data0_av, displ2_data1_av, displ2_data2_av; //0=X, 1=Y, 2=Z double displ2_data0_var, displ2_data1_var, displ2_data2_var; unsigned long long count_num_datapoints = 0; { //calc displ2 bucket avg //unsigned long long displ2_data0, displ2_data1, displ2_data2; //temp counters //0=X, 1=Y, 2=Z displ2_data0_av = displ2_data1_av = displ2_data2_av = 0; for (i=0; i<total_num_buckets_displ2; i++) { //displ2_data0_av += (i+.5)*total_num_buckets_displ2*ds_cube_dist_displ2->data[0][bucket_index][i]; //displ2_data1_av += ds_cube_dist_displ2->data[1][bucket_index][i]; //displ2_data2_av += ds_cube_dist_displ2->data[2][bucket_index][i]; displ2_data0_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[0][bucket_index][i]; displ2_data1_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[1][bucket_index][i]; displ2_data2_av += ((double)i+.5) * (double)ds_cube_dist_displ2->data[2][bucket_index][i]; count_num_datapoints += ds_cube_dist_displ2->data[0][bucket_index][i]; //sum should be equal for x, y and z data } //displ2_data0_av /= (double)displ2_data0_av / (double)total_num_buckets_displ2; //displ2_data1_av /= (double)displ2_data1_av / (double)total_num_buckets_displ2; //displ2_data2_av /= (double)displ2_data2_av / (double)total_num_buckets_displ2; displ2_data0_av /= (double)count_num_datapoints; displ2_data1_av /= (double)count_num_datapoints; displ2_data2_av /= (double)count_num_datapoints; } //rw_raw_data_fout << "displ2 bucket av_, _x, _y, _z" << endl; //rw_raw_data_fout << ", " << displ2_data0_av; //rw_raw_data_fout << ", " << displ2_data1_av; //rw_raw_data_fout << ", " << displ2_data2_av; //rw_raw_data_fout << endl; ////calc displ2 bucket standard dev //displ2_data0_sd = displ2_data1_sd = displ2_data2_sd = 0; //for (i=0; i<total_num_buckets_displ2; i++) //{ // displ2_data0_sd += ds_cube_dist_displ2->data[0][bucket_index][i]; // displ2_data1_sd += ds_cube_dist_displ2->data[1][bucket_index][i]; // displ2_data2_sd += ds_cube_dist_displ2->data[2][bucket_index][i]; //} //calc displ2 bucket variance displ2_data0_var = displ2_data1_var = displ2_data2_var = 0; for (i=0; i<total_num_buckets_displ2; i++) { displ2_data0_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[0][bucket_index][i] - displ2_data0_av, 2); displ2_data1_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[1][bucket_index][i] - displ2_data1_av, 2); displ2_data2_var += pow( (double)(i+.5) * (double)ds_cube_dist_displ2->data[2][bucket_index][i] - displ2_data2_av, 2); } displ2_data0_av *= buckets_displ2_width; displ2_data1_av *= buckets_displ2_width; displ2_data2_av *= buckets_displ2_width; displ2_data0_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; displ2_data1_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; displ2_data2_var *= pow(buckets_displ2_width, 2) / (double)count_num_datapoints; rw_raw_data_fout << "displ2 av_,"; rw_raw_data_fout << ", " << displ2_data0_av; rw_raw_data_fout << ", " << displ2_data1_av; rw_raw_data_fout << ", " << displ2_data2_av; rw_raw_data_fout << endl; rw_raw_data_fout << "displ2 sd_,"; rw_raw_data_fout << ", " << sqrt(displ2_data0_var); rw_raw_data_fout << ", " << sqrt(displ2_data1_var); rw_raw_data_fout << ", " << sqrt(displ2_data2_var); rw_raw_data_fout << endl; rw_raw_data_fout << "displ2 var_,"; rw_raw_data_fout << ", " << displ2_data0_var; rw_raw_data_fout << ", " << displ2_data1_var; rw_raw_data_fout << ", " << displ2_data2_var; rw_raw_data_fout << endl; //estimate gamma distrib parameters double alphax, betax, alphay, betay, alphaz, betaz; alphax = pow(displ2_data0_av, 2) / displ2_data0_var; betax = displ2_data0_var / displ2_data0_av; alphay = pow(displ2_data1_av, 2) / displ2_data1_var; betay = displ2_data1_var / displ2_data1_av; alphaz = pow(displ2_data2_av, 2) / displ2_data2_var; betaz = displ2_data2_var / displ2_data2_av; rw_raw_data_fout << ", k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; rw_raw_data_fout << endl; rw_raw_data_fout << ", " << alphax << ", " << betax; rw_raw_data_fout << ", " << alphay << ", " << betay; rw_raw_data_fout << ", " << alphaz << ", " << betaz; rw_raw_data_fout << endl; ////save whole histogram data to file //for (i=0; i<total_num_buckets_displ2; i++) //{ // rw_raw_data_fout << i << ","; // rw_raw_data_fout << ds_cube_dist_displ2->data[0][bucket_index][i] << ","; // rw_raw_data_fout << ds_cube_dist_displ2->data[1][bucket_index][i] << ","; // rw_raw_data_fout << ds_cube_dist_displ2->data[2][bucket_index][i] << ","; // rw_raw_data_fout << endl; //} unsigned long long total_ds_cube_dist_displ2_data0, total_ds_cube_dist_displ2_data1, total_ds_cube_dist_displ2_data2; total_ds_cube_dist_displ2_data0 = 0; total_ds_cube_dist_displ2_data1 = 0; total_ds_cube_dist_displ2_data2 = 0; //save whole histogram data to file + add Gamma distribution values for (i=0; i<total_num_buckets_displ2; i++) { double t = (double)(i+.5)*(double)buckets_displ2_width; rw_raw_data_fout << i << ","; rw_raw_data_fout << t << ","; //middle of displ2 bucket av //rw_raw_data_fout << << ","; //true av for displ2_x is not available here - would require extra structure to hold data rw_raw_data_fout << ds_cube_dist_displ2->data[0][bucket_index][i] << ","; total_ds_cube_dist_displ2_data0 += ds_cube_dist_displ2->data[0][bucket_index][i]; rw_raw_data_fout << ds_cube_dist_displ2->data[1][bucket_index][i] << ","; total_ds_cube_dist_displ2_data1 += ds_cube_dist_displ2->data[1][bucket_index][i]; rw_raw_data_fout << ds_cube_dist_displ2->data[2][bucket_index][i] << ","; total_ds_cube_dist_displ2_data2 += ds_cube_dist_displ2->data[2][bucket_index][i]; rw_raw_data_fout << gamma_pdf(t, alphax, betax) << ","; rw_raw_data_fout << gamma_pdf(t, alphay, betay) << ","; rw_raw_data_fout << gamma_pdf(t, alphaz, betaz) << ","; rw_raw_data_fout << endl; } rw_raw_data_fout << "total,,"; rw_raw_data_fout << total_ds_cube_dist_displ2_data0 << ","; rw_raw_data_fout << total_ds_cube_dist_displ2_data1 << ","; rw_raw_data_fout << total_ds_cube_dist_displ2_data2 << ","; rw_raw_data_fout << endl; //file sub-footer / per dist bucket rw_raw_data_fout << "*****************************************************************************"; rw_raw_data_fout << endl << endl; } //close raw rw file rw_raw_data_fout << endl;// << "****" << endl << endl; rw_raw_data_fout.close(); //delete histogram data cube delete ds_cube_dist_displ2; delete ds_cube_min_displ2_per_dist_bucket; } else {//dont save full data while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); count_distance_bucket[bucket_index] ++; distance_bucket[bucket_index] += (double)temp_entry.distance; disp2x_bucket[bucket_index] += (double)temp_entry.dispx; disp2y_bucket[bucket_index] += (double)temp_entry.dispy; disp2z_bucket[bucket_index] += (double)temp_entry.dispz; } //no save_full_RW_data, next datapoint in the distance structure distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } } // should depend on total_num_buckets? // depending on a fixed value -> not good //min_required_num_per_bucket = (unsigned int)(count_distances*0.02); //xx% of total # of datapoints -> not good either // depending on value in most populated bucket for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) if ( max_num_in_a_bucket < count_distance_bucket[bucket_index] ) max_num_in_a_bucket = count_distance_bucket[bucket_index]; min_required_num_per_bucket = (unsigned int)(max_num_in_a_bucket*percent_densest_bucket); //percent_densest_bucket% of most populated bucket num_used_buckets = 0; count_distance_bucket_for_grand_av = 0; total_count_distance_bucket = 0; SSxx_distance = SSxx_distance_for_grand_av = (double)0; SSyy_dispx = SSyy_dispy = SSyy_dispz = (double)0; SSxy_ddx = SSxy_ddy = SSxy_ddz = (double)0; for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) {//summing everything - SS is just the overall SUM here -> temp variables to calc overall AVERAGES first if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) //{//keep value only if there are enough datapoints in that particular bucket // num_used_buckets ++; // SSxx_distance += distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; // total_count_distance_bucket += count_distance_bucket[bucket_index]; // //SSxx_distance += distance_bucket[bucket_index] = ((double)bucket_index+.5) * (double)bucket_width; // //SSxx_distance calc ???should??? involve middle of bucket value, instead of real data value -> done // SSyy_dispx += disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; // SSyy_dispy += disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; // SSyy_dispz += disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //} {//keep value only if there are enough datapoints in that particular bucket //this is to calculate the grand_av - expected to be much smaller than the lin reg av SSxx_distance_for_grand_av += distance_bucket[bucket_index]; SSxy_ddx += disp2x_bucket[bucket_index]; SSxy_ddy += disp2y_bucket[bucket_index]; //these are not SSxys they are temp variables reused to calc the grand_av SSxy_ddz += disp2z_bucket[bucket_index]; count_distance_bucket_for_grand_av += count_distance_bucket[bucket_index]; //SSxx_distance += distance_bucket[bucket_index] = ((double)bucket_index+.5) * (double)bucket_width; //SSxx_distance calculation shouldn't involve middle of bucket value, instead better to use real data value: SSxx_distance += distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //this is calculating avrg for lin regression of ensemble averages (ie per dist bucket - not grand_av) SSyy_dispx += disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispy += disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; SSyy_dispz += disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //av per bucket - calc'd above already //distance_bucket[bucket_index] = distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2x_bucket[bucket_index] = disp2x_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2y_bucket[bucket_index] = disp2y_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; //disp2z_bucket[bucket_index] = disp2z_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]; num_used_buckets ++; total_count_distance_bucket += count_distance_bucket[bucket_index]; } else { distance_bucket[bucket_index] = 0; disp2x_bucket[bucket_index] = 0; //averages per bucket disp2y_bucket[bucket_index] = 0; disp2z_bucket[bucket_index] = 0; } } percent_densest_bucket -= (float).01; //decrease % to ensure at least 3 buckets } while ( (num_used_buckets < 3) && (percent_densest_bucket >= 0) ); if (num_used_buckets == 1 ) { write_to_log("only 1 bucket!"); } ////calc overall averages //av_distance = (float)SSxx_distance/(float)num_used_buckets; //xbar //av_dispx = (float)SSyy_dispx/(float)num_used_buckets; //av_dispy = (float)SSyy_dispy/(float)num_used_buckets; //ybar //av_dispz = (float)SSyy_dispz/(float)num_used_buckets; //calc grand (overall) averages av_distance = (float)(SSxx_distance_for_grand_av/(double)count_distance_bucket_for_grand_av); //grand xbar av_dispx = (float)(SSxy_ddx/(double)count_distance_bucket_for_grand_av); av_dispy = (float)(SSxy_ddy/(double)count_distance_bucket_for_grand_av); //grand ybar av_dispz = (float)(SSxy_ddz/(double)count_distance_bucket_for_grand_av); stats_results_fout << "grand avg, " << endl; stats_results_fout << "dist, displ2X, displ2Y, displ2Z, " << endl; stats_results_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; //estimated D/D0 coefficients stats_results_fout << "and estimated D/D0 coefficients (xyz)" << endl; { D_D0x = (float)1.5*av_dispx/av_distance/lambda_micron; D_D0y = (float)1.5*av_dispy/av_distance/lambda_micron; D_D0z = (float)1.5*av_dispz/av_distance/lambda_micron; } stats_results_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << endl; stats_results_fout << endl; //calc lin reg avg av_distance = (float)SSxx_distance/(float)num_used_buckets; //xbar av_dispx = (float)SSyy_dispx/(float)num_used_buckets; av_dispy = (float)SSyy_dispy/(float)num_used_buckets; //ybar av_dispz = (float)SSyy_dispz/(float)num_used_buckets; stats_results_fout << "lin reg avg, " << endl; stats_results_fout << "dist, displ2X, displ2Y, displ2Z, " << endl; stats_results_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; //estimated D/D0 coefficients stats_results_fout << "and estimated D/D0 coefficients (xyz)" << endl; { D_D0x = (float)1.5*av_dispx/av_distance/lambda_micron; D_D0y = (float)1.5*av_dispy/av_distance/lambda_micron; D_D0z = (float)1.5*av_dispz/av_distance/lambda_micron; } stats_results_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << endl; stats_results_fout << endl; //now that we have the averages //calc standard errors for everything // Grab start of list pointer distance_structure_list_node_ptr = distance_structure_list_ptr->getHeadPtr(); while ( distance_structure_list_node_ptr ) {//summing everything - _bucket is just the SUM here -> to calc AVERAGES first distance_structure temp_entry; temp_entry = distance_structure_list_node_ptr->item; if (temp_entry.distance < limited_max_distance) { bucket_index = (unsigned int)(temp_entry.distance/bucket_width); if (count_distance_bucket[bucket_index] > min_required_num_per_bucket) { sd_distance_bucket[bucket_index] += pow((double)temp_entry.distance - distance_bucket[bucket_index], 2); sd_dispx_bucket[bucket_index] += pow((double)temp_entry.dispx - disp2x_bucket[bucket_index], 2); sd_dispy_bucket[bucket_index] += pow((double)temp_entry.dispy - disp2y_bucket[bucket_index], 2); sd_dispz_bucket[bucket_index] += pow((double)temp_entry.dispz - disp2z_bucket[bucket_index], 2); } } distance_structure_list_node_ptr = distance_structure_list_node_ptr->next; } //sub part header stats_results_fout << "(dist in microns), (displ^2 in microns^2)" << endl; stats_results_fout << "av_dist, av_displ2X, av_displ2Y, av_displ2Z, count, "; stats_results_fout << "sd_dist_norm'd, sd_disp2X_norm'd, sd_disp2Y_norm'd, sd_disp2Z_norm'd, "; //stats_results_per_file_number_fout << "(av-sd)_dist, (av+sd)_dist, (av-sd)_disp2X, (av+sd)_disp2X, "; //stats_results_per_file_number_fout << "(av-sd)_disp2Y, (av+sd)_disp2Y, (av-sd)_disp2Z, (av+sd)_disp2Z, "; //stats_results_per_file_number_fout << "k=alphaX, theta=betaX, k=alphaY, theta=betaY, k=alphaZ, theta=betaZ, "; stats_results_fout << endl; //init variables SSxx_distance = (double)0; //SSxx SSyy_dispx = (double)0; SSyy_dispy = (double)0; //SSyy SSyy_dispz = (double)0; SSxy_ddx = (double)0; SSxy_ddy = (double)0; //SSxy SSxy_ddz = (double)0; //calc SS for (bucket_index=0; bucket_index<total_num_buckets; bucket_index++) if (distance_bucket[bucket_index] > 0) { SSxx_distance += pow(distance_bucket[bucket_index], 2); //SSxx=SUM(x^2) SSyy_dispx += pow(disp2x_bucket[bucket_index], 2); SSyy_dispy += pow(disp2y_bucket[bucket_index], 2); //SSyy=SUM(y^2) SSyy_dispz += pow(disp2z_bucket[bucket_index], 2); SSxy_ddx += distance_bucket[bucket_index] * disp2x_bucket[bucket_index]; SSxy_ddy += distance_bucket[bucket_index] * disp2y_bucket[bucket_index]; //SSxy=SUM(x*y) SSxy_ddz += distance_bucket[bucket_index] * disp2z_bucket[bucket_index]; sd_distance_bucket[bucket_index] = sqrt(sd_distance_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispx_bucket[bucket_index] = sqrt(sd_dispx_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispy_bucket[bucket_index] = sqrt(sd_dispy_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); sd_dispz_bucket[bucket_index] = sqrt(sd_dispz_bucket[bucket_index]/(double)count_distance_bucket[bucket_index]); //write values to file - to check the lin. regression fit stats_results_fout << distance_bucket[bucket_index] << ", " \ << disp2x_bucket[bucket_index] << ", " \ << disp2y_bucket[bucket_index] << ", " \ << disp2z_bucket[bucket_index] << ", " \ << count_distance_bucket[bucket_index] << ", " \ << (sd_distance_bucket[bucket_index] / distance_bucket[bucket_index]) << ", " \ << (sd_dispx_bucket[bucket_index] / disp2x_bucket[bucket_index]) << ", " \ << (sd_dispy_bucket[bucket_index] / disp2y_bucket[bucket_index]) << ", " \ << (sd_dispz_bucket[bucket_index] / disp2z_bucket[bucket_index]) << ", "; //<< (distance_bucket[bucket_index] - sd_distance_bucket[bucket_index]) << ", " \ //<< (distance_bucket[bucket_index] + sd_distance_bucket[bucket_index]) << ", " \ //<< (disp2x_bucket[bucket_index] - sd_dispx_bucket[bucket_index]) << ", " \ //<< (disp2x_bucket[bucket_index] + sd_dispx_bucket[bucket_index]) << ", " \ //<< (disp2y_bucket[bucket_index] - sd_dispy_bucket[bucket_index]) << ", " \ //<< (disp2y_bucket[bucket_index] + sd_dispy_bucket[bucket_index]) << ", " \ //<< (disp2z_bucket[bucket_index] - sd_dispz_bucket[bucket_index]) << ", " \ //<< (disp2z_bucket[bucket_index] + sd_dispz_bucket[bucket_index]) << ", "; stats_results_fout << endl; } delete [] count_distance_bucket; delete [] distance_bucket; delete [] disp2x_bucket; delete [] disp2y_bucket; delete [] disp2z_bucket; delete [] sd_distance_bucket; delete [] sd_dispx_bucket; delete [] sd_dispy_bucket; delete [] sd_dispz_bucket; //save results for 0-intercept float Sumxx, Sumxyx, Sumxyy, Sumxyz; Sumxx = (float)SSxx_distance; Sumxyx = (float)SSxy_ddx; Sumxyy = (float)SSxy_ddy; Sumxyz = (float)SSxy_ddz; //finally: SSxx_distance = SSxx_distance - (double)num_used_buckets * pow((double)av_distance, 2); //SSxx=SUM(x^2)-n*xbar^2 SSyy_dispx = SSyy_dispx - (double)num_used_buckets * pow((double)av_dispx, 2); SSyy_dispy = SSyy_dispy - (double)num_used_buckets * pow((double)av_dispy, 2); //SSyy=SUM(y^2)-n*ybar^2 SSyy_dispz = SSyy_dispz - (double)num_used_buckets * pow((double)av_dispz, 2); SSxy_ddx = SSxy_ddx - (double)num_used_buckets * (double)av_distance * (double)av_dispx; SSxy_ddy = SSxy_ddy - (double)num_used_buckets * (double)av_distance * (double)av_dispy; //SSxy=SUM(x*y)-n*xbar*ybar SSxy_ddz = SSxy_ddz - (double)num_used_buckets * (double)av_distance * (double)av_dispz; //get a, b (y=a+bx) -- disp^2 vs distance //b=SSxy/SSxx float bx, by, bz; bx = (float)(SSxy_ddx/SSxx_distance); by = (float)(SSxy_ddy/SSxx_distance); bz = (float)(SSxy_ddz/SSxx_distance); //(all we want is the slope b) //now a=ybar-b*xbar float ax, ay, az; ax = av_dispx - bx*av_distance; ay = av_dispy - by*av_distance; az = av_dispz - bz*av_distance; float sigmax_distance, sigmay_dispx, sigmay_dispy, sigmay_dispz, covxy_ddx, covxy_ddy, covxy_ddz; //sigmax^2=SSxx/n sigmax_distance = (float)(SSxx_distance / num_used_buckets); //sigmay^2=SSyy/n sigmay_dispx = (float)(SSyy_dispx / (double)num_used_buckets); sigmay_dispy = (float)(SSyy_dispy / (double)num_used_buckets); sigmay_dispz = (float)(SSyy_dispz / (double)num_used_buckets); //Cov(x,y)=SSxy/n covxy_ddx = (float)(SSxy_ddx / (double)num_used_buckets); covxy_ddy = (float)(SSxy_ddy / (double)num_used_buckets); covxy_ddz = (float)(SSxy_ddz / (double)num_used_buckets); //s=sqrt( (SSyy-b*SSxy)/(n-2) ) float s_x, s_y, s_z; s_x = (float)sqrt( (SSyy_dispx-bx*SSxy_ddx)/((float)num_used_buckets-2) ); s_y = (float)sqrt( (SSyy_dispy-by*SSxy_ddy)/((float)num_used_buckets-2) ); s_z = (float)sqrt( (SSyy_dispz-bz*SSxy_ddz)/((float)num_used_buckets-2) ); //se(a)=s*sqrt( 1/n + xbar^2/SSxx ) float se_ax, se_ay, se_az; se_ax = s_x*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_ay = s_y*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); se_az = s_z*(float)sqrt( 1/(float)num_used_buckets + av_distance*av_distance/SSxx_distance); //se(b)=s/sqrt(SSxx) float se_bx, se_by, se_bz; se_bx = s_x / (float)sqrt(SSxx_distance); se_by = s_y / (float)sqrt(SSxx_distance); se_bz = s_z / (float)sqrt(SSxx_distance); //R^2=SSxy^2/SSxx/SSyy float R2_x, R2_y, R2_z; R2_x = (float)(SSxy_ddx*SSxy_ddx/SSxx_distance/SSyy_dispx); R2_y = (float)(SSxy_ddy*SSxy_ddy/SSxx_distance/SSyy_dispy); R2_z = (float)(SSxy_ddz*SSxy_ddz/SSxx_distance/SSyy_dispz); //store stats results to file stats_results_fout << "average" << endl; stats_results_fout << av_distance << ", " << av_dispx << ", " << av_dispy << ", " << av_dispz << endl; stats_results_fout << "sigma^2" << endl; stats_results_fout << sigmax_distance << ", " << sigmay_dispx << ", " << sigmay_dispy << ", " << sigmay_dispz << endl << endl; stats_results_fout << "(y=a+bx) -- disp^2 vs distance" << endl; stats_results_fout << ", disp2x, disp2y, disp2z, " << endl; stats_results_fout << "a, " << ax << ", " << ay << ", " << az << endl; stats_results_fout << "b, " << bx << ", " << by << ", " << bz << endl; stats_results_fout << "R^2, " << R2_x << ", " << R2_y << ", " << R2_z << endl; stats_results_fout << "se_a, " << se_ax << ", " << se_ay << ", " << se_az << endl; stats_results_fout << "se_b, " << se_bx << ", " << se_by << ", " << se_bz << endl; //null intercept, b=Sumxy[]/Sumxx stats_results_fout << "(y=bx) -- zero (0) intercept" << endl; stats_results_fout << "b, " << (float)(Sumxyx/Sumxx) \ << ", " << (float)(Sumxyy/Sumxx) \ << ", " << (float)(Sumxyz/Sumxx) << endl; //now get D/D0 stats_results_fout << endl; //stats_results_fout << "D/D0 coefficients (xyz),,,,,coeff*porosity" << endl; stats_results_fout << "D/D0 coefficients (xyz)," << endl; { D_D0x = (float)1.5*bx/lambda_micron; D_D0y = (float)1.5*by/lambda_micron; D_D0z = (float)1.5*bz/lambda_micron; } stats_results_fout << ", " << D_D0x << ", " << D_D0y << ", " << D_D0z << ", " << endl; //stats_results_fout << ", " << D_D0x*average_ROI_porosity << ", " << D_D0y*average_ROI_porosity << ", " << D_D0z*average_ROI_porosity << endl; stats_results_fout << "+/-" << endl; stats_results_fout << ", " << 1.5*se_bx/lambda_micron \ << ", " << 1.5*se_by/lambda_micron \ << ", " << 1.5*se_bz/lambda_micron << endl; stats_results_fout << "with zero (0) intercept" << endl; stats_results_fout << ", " << 1.5*(float)(Sumxyx/Sumxx)/lambda_micron \ << ", " << 1.5*(float)(Sumxyy/Sumxx)/lambda_micron \ << ", " << 1.5*(float)(Sumxyz/Sumxx)/lambda_micron << endl; stats_results_fout << endl; //compare to known literature values stats_results_fout << "D/D0 coefficients from literature" << endl; stats_results_fout << "porosity, Maxwell, Wakao, Rayleigh, " << endl; stats_results_fout << average_ROI_porosity << ", " ; stats_results_fout << 2*average_ROI_porosity/(3-average_ROI_porosity) << ", " ; stats_results_fout << average_ROI_porosity*average_ROI_porosity<< ", " ; // /sqrt((float)3) << ", " ; float rayleigh_result; rayleigh_result = (float).3938*pow((float)1-average_ROI_porosity, (float)10/(float)3); stats_results_fout << (2*average_ROI_porosity-rayleigh_result)/(3-average_ROI_porosity-rayleigh_result) << ", " ; } else //(rod==0) { if (!rod_zero_aborted_enough_data_for_stats) //(dlg->stop_request) { //if (total_count_distance_bucket < 500000) write_to_log("no results available: not enough data"); goto delvar; } else {//enough data, can proceed float average_ROI_porosity = total_ROI_porosity/(float)total_file_number; stats_results_fout << "av_true_ROI_porosity, " << average_ROI_porosity << endl; stats_results_fout << "av_ROI_porosity_from tracer start, " \ << (float)total_count_tracer_start_hit_pore/(float)total_count_tracer_start_hit << endl; coef_rod_x /= (float)num_coeff; coef_rod_y /= (float)num_coeff; coef_rod_z /= (float)num_coeff; coef_rod_sex /= (float)num_coeff; coef_rod_sey /= (float)num_coeff; coef_rod_sez /= (float)num_coeff; coef_rod0_x /= (float)num_coeff; coef_rod0_y /= (float)num_coeff; coef_rod0_z /= (float)num_coeff; stats_results_fout << "average_D/D0 coefficients (xyz)" << endl; stats_results_fout << coef_rod_x << ", " ; stats_results_fout << coef_rod_y << ", " ; stats_results_fout << coef_rod_z << ", " << endl ; stats_results_fout << "average_SE of D/D0 coefficients (xyz)" << endl; stats_results_fout << coef_rod_sex << ", " ; stats_results_fout << coef_rod_sey << ", " ; stats_results_fout << coef_rod_sez << ", " << endl ; stats_results_fout << "average_D/D0 coefficients - 0 intercept(xyz)" << endl; stats_results_fout << coef_rod0_x << ", " ; stats_results_fout << coef_rod0_y << ", " ; stats_results_fout << coef_rod0_z << ", " << endl ; stats_results_fout << endl; stats_results_fout << "D/D0 coefficients from literature" << endl; stats_results_fout << "porosity, Maxwell, Wakao, Rayleigh, " << endl; stats_results_fout << average_ROI_porosity << ", " ; stats_results_fout << 2*average_ROI_porosity/(3-average_ROI_porosity) << ", " ; stats_results_fout << average_ROI_porosity*average_ROI_porosity << ", " ;// /sqrt((float)3) << ", " ; float rayleigh_result; rayleigh_result = (float).3938*pow((float)1-average_ROI_porosity, (float)10/(float)3); stats_results_fout << (2*average_ROI_porosity-rayleigh_result)/(3-average_ROI_porosity-rayleigh_result) << ", " ; } } //close stats values file stats_results_fout << endl << "****" << endl << endl; stats_results_fout.close(); delvar: delete distMap; if (!dlg->stop_request) { write_to_log("Simultaneous Diffusion - done"); dlg->txt_field_progress.SetWindowText("Simultaneous Diffusion - done"); } else { //close stats values file if (stats_results_fout.is_open()) { stats_results_fout << endl << "ABORTED" << endl << endl; stats_results_fout.close(); } //close file for saving full RW raw data if (rw_raw_data_fout.is_open()) { rw_raw_data_fout << endl << "ABORTED" << endl << endl; rw_raw_data_fout.close(); } //close file for saving cv data if (check_cv) { rw_cv_data_fout << endl << "ABORTED @ total_particle_end_of_run, " << total_particle_end_of_run << endl << endl; rw_cv_data_fout.close(); } write_to_log("Simultaneous Diffusion - ABORTED"); dlg->txt_field_progress.SetWindowText("Simultaneous Diffusion - ABORTED"); //should also dist_n_disp_sq_fout.close();??? //if (dist_n_disp_sq_fout.is_open()) // dist_n_disp_sq_fout.close(); } delete distance_structure_list_ptr; }