/* Brief explanation: cDSinfo* ds_cube contains a binarized porous structure cube (compressed... 1 bit for each voxel). Then 'nPores' are run through this 'ds' structure in the following manner. First, a random porous spot on the top layer is picked as the starting point. Then the tracer tries going straight down. If the path is blocked, it tries all 8 down+left/right (diagonal) combinations If this fails, it tries going into one of the 8 sideways directions Finally, it will try to go back up/diagonally The pore cannot go back through the path it travelled before (an infinite loop will occur if this is allowed), so if it gets trapped, that particular run is rejected and a new one is begun. */ #include "stdafx.h" #include "Tortuosity.h" using namespace std; /*---------------- GENERATES THE RANDOM NUMBER -----------------*/ int randomNo(int k) {//should be FLOAT return(1 + rand()*k/(RAND_MAX + 1)); } /* TORTUOSITY MEMBERS */ Tortuosity::Tortuosity(DS_1b_cube *ds_new_cube, void* new_TortDlg) { ds_cube = ds_new_cube; wZ = ds_cube->wZ; wY = ds_cube->wY; wX = ds_cube->wX; TortDlg = (TortuosityDialog*) new_TortDlg; path_cube = new DS_1b_cube(wZ, wY, wX); head = new vector_node; head->next = NULL; tail = head; } Tortuosity::~Tortuosity() { delete path_cube; delete head; } void Tortuosity::add_vector_node(unsigned short int z, unsigned short int y, unsigned short int x) { vector_node * new_node = new vector_node;//(vector_node*) malloc(sizeof(vector_node)); new_node->next = NULL; new_node->x = x; new_node->y = y; new_node->z = z; tail->next = new_node; tail = new_node; } void Tortuosity::clear_vector_list() { tail = head->next; delete head; //erase linked list one spot at a time while(tail) { path_cube->data[tail->z][tail->y][(tail->x)/8] = 0x00; head = tail->next; delete tail; tail = head; } head = new vector_node; head->next = NULL; tail = head; } /*----------- MAIN FUNCTION TO CALL -------------*/ /* bmp_name, sX, sY, res, manual, invert, T are only used for logs startingSlice, lastSlice used to calc N, and logs return values: 0 successful 1 failed, could not find a valid start location after 100000 tries 2 failed, could not find a complete path through structure after 100000 tries 3 file opening error 4 memory allocation problem */ unsigned short int Tortuosity::run_tortuosity(char* out_name, int nPores, unsigned short int depth, int rec_tracer_num){ //global variables double update_freq = nPores/100; double curr_bar = 0; double tracer_cpy_freq=0; double curr_cpy_cnt=0; int poresFound; double Tortuosity; unsigned short int i, j, k; int failed_attempts; char file_name[512]; FILE *fout1, *tracer_file; /* if user wants tracers saved, also create a cube for that*/ if (rec_tracer_num > 0){ tracer_cpy_freq = (nPores)/(rec_tracer_num); curr_cpy_cnt = tracer_cpy_freq - 1; tracers_cube = new DS_1b_cube(wZ, wY, wX); if (!tracers_cube) { return 4; } } // clear path cube path_cube->clear(); /* create a data output file */ sprintf(file_name,"%s_data.txt",out_name); fout1 = fopen(file_name, "w"); if ( fout1 == NULL ) { return 3; } /*--------------- CODE STARTS AT 'X' EQUAL TO ZERO ----------*/ poresFound=0; failed_attempts=0; //this will run one tracer down every open pore if (nPores < 1){ //find number of valid spots for (j = 0; j < wY; j++){ for (k = 0; k < wX; k++){ if (!ds_cube->get_spot(0, j, k) ) poresFound++; } } TortDlg->SetIntValField(&(TortDlg->field_run_num_XY_cntrl), poresFound); update_freq = poresFound/100; if (rec_tracer_num > 0){ tracer_cpy_freq = poresFound/(rec_tracer_num); curr_cpy_cnt = tracer_cpy_freq - 1; } poresFound = 0; for (j = 0; j < wY; j++){ for (k = 0; k < wX; k++){ if (!ds_cube->get_spot(0, j, k) ){ poresFound++; length = 0; s = 0; r = j; c = k; //mark the starting spot as traveled path_cube->set_spot_1(0, j, k); if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt) tracers_cube->set_spot_1(0, j, k); /* travel through structure until last slice reached */ while (s < (depth-1) ){ //if the tracer got stuck in a dead end, restart process if (!send_tracer()){ failed_attempts++; goto skipToTheEnd2; } // if user requested tracers to be saved, update the new voxel if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt) tracers_cube->set_spot_1(s, r, c); } /* update recording of tracers */ if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt){ curr_cpy_cnt = curr_cpy_cnt + tracer_cpy_freq; } /* update GUI progress bar */ if (poresFound > curr_bar){ curr_bar = curr_bar + update_freq; TortDlg->pbar_simulation_progress.StepIt(); } //compute tortuosity Tortuosity = length/(wZ-1); //store to result file... "tracer_number tortuosity" fprintf(fout1,"%d \t %f\n", poresFound, Tortuosity); nPores++; //number of successful pores skipToTheEnd2: ; } } } TortDlg->SetIntValField(&(TortDlg->field_run_num_XY_cntrl), nPores); } //nPores created in random start locations else{ while (poresFound < nPores){ s=0; length=0; //find a starting location (repeat until an an empty pore is found) do { /* generate a 0 <= random < width */ r = (unsigned short)randomNo(wY-1); c = (unsigned short)randomNo(wX-1); /* in case of a bad input, don't run forever */ failed_attempts++; if (failed_attempts > 100000) { /* clean up */ if (rec_tracer_num > 0) delete tracers_cube; fclose(fout1); return 1; } } while(ds_cube->get_spot(s, r, c) ); //mark the starting spot as traveled path_cube->set_spot_1(s, r, c); if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt) tracers_cube->set_spot_1(s, r, c); /* travel through structure until last slice reached */ while (s < (depth-1) ){ //add new node to list of traveled locations add_vector_node(s, r, c); //if the tracer got stuck in a dead end, restart process if (!send_tracer()){ failed_attempts++; goto skipToTheEnd; } // if user requested tracers to be saved, update the new voxel if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt) tracers_cube->set_spot_1(s, r, c); } path_cube->data[s][r][(c)/8] = 0x00; /* update recording of tracers */ if (rec_tracer_num > 0) if (poresFound >= curr_cpy_cnt){ curr_cpy_cnt = curr_cpy_cnt + tracer_cpy_freq; } /* update GUI progress bar */ if (poresFound > curr_bar){ curr_bar = curr_bar + update_freq; TortDlg->pbar_simulation_progress.StepIt(); } //increase number of successful tracers, reset attempts poresFound++; failed_attempts = 0; //compute tortuosity Tortuosity = length/(wZ-1); //store to result file... "tracer_number tortuosity" fprintf(fout1,"%d \t %f\n", poresFound, Tortuosity); skipToTheEnd: clear_vector_list(); /* fail if 100000 or more unsuccessful tries */ if (failed_attempts > 100000) { /* clean up */ if (rec_tracer_num > 0) delete tracers_cube; fclose(fout1); return 2; } } } //fprintf(foutlog, "Run successfully completed! See file %s for individual run results.\n", results_file); fflush(foutlog); write_to_log("Tortuosity: All tracers have been successfully sent through structure."); if (rec_tracer_num > 0){ sprintf(file_name, "%s_saved_tracers.stack.bin", out_name); tracer_file = fopen(file_name,"w"); for (k=0; k < depth; k++){ for (i = 0; i < wY; i++){ for (j = 0; j < wX; j++){ if (!tracers_cube->get_spot(k, i, j)) fputc(0x00,tracer_file); else fputc(0xFF,tracer_file); } } } fclose(tracer_file); delete tracers_cube; } //clean up fclose(fout1); return 0; } /*------------------ TRACE FUNCTION BEGINS HERE ---------------*/ /* First the tracer tries going straight down. If the path is blocked, it tries all 8 down+left/right (diagonal) combinations If this fails, it tries going into one of the 8 sideways directions Finally, it will try to go back up/diagonally/* */ bool Tortuosity::send_tracer() { /*---------------- FIRST PHASE OF LOWER SLICE ---------------*/ if (ds_cube->get_spot((s+1), r, c) == false) { if(path_cube->set_spot_1((s+1), (r), (c))==true){ length=length+1; s=s+1; return true; } } /*--------------- SECOND PHASE OF LOWER SLICE ---------------*/ if ( (r-1) >= 0 ) { if (ds_cube->get_spot((s+1), (r-1), c) == false) { if(path_cube->set_spot_1((s+1), (r-1), c)==true){ length=length+1.414; s=s+1; r=r-1; return true; } } } if ( (c-1) >= 0 ) { if (ds_cube->get_spot((s+1), r, (c-1)) == false) { if(path_cube->set_spot_1((s+1), r, (c-1))==true){ length=length+1.414; s=s+1; c=c-1; return true; } } } if ( (c+1) < wX ) { if (ds_cube->get_spot((s+1), r, (c+1)) == false) { if(path_cube->set_spot_1((s+1), r, (c+1))==true){ length=length+1.414; s=s+1; c=c+1; return true; } } } if ( (r+1) < wY ) { if (ds_cube->get_spot((s+1), (r+1), c) == false) { if(path_cube->set_spot_1((s+1), (r+1), c)==true){ length=length+1.414; s=s+1; r=r+1; return true; } } } /*--------------- THIRD PHASE OF LOWER SLICE ---------------*/ if ( ((r-1) >= 0) && ( (c-1) >= 0) ){ if (ds_cube->get_spot((s+1), (r-1), (c-1)) == false) { if(path_cube->set_spot_1((s+1), (r-1), (c-1))==true){ length=length+1.732; s=s+1; r=r-1; c=c-1; return true; } } } if ( ((r-1) >= 0) && ( (c+1) < wX) ){ if (ds_cube->get_spot((s+1), (r-1), (c+1)) == false) { if(path_cube->set_spot_1((s+1), (r-1), (c+1))==true){ length=length+1.732; s=s+1; r=r-1; c=c+1; return true; } } } if ( ((r+1) < wY) && ( (c-1) >= 0) ){ if (ds_cube->get_spot((s+1), (r+1), (c-1)) == false) { if(path_cube->set_spot_1((s+1), (r+1), (c-1))==true){ length=length+1.732; s=s+1; r=r+1; c=c-1; return true; } } } if ( ((r+1) < wY) && ( (c+1) < wX) ){ if (ds_cube->get_spot((s+1), (r+1), (c+1)) == false) { if(path_cube->set_spot_1((s+1), (r+1), (c+1))==true){ length=length+1.732; s=s+1; r=r+1; c=c+1; return true; } } } /*----------- FIRST PHASE OF SAME SLICE -----------------*/ if ( (r-1) >= 0){ if (ds_cube->get_spot(s, (r-1), c) == false) { if(path_cube->set_spot_1((s), (r-1), (c))==true){ length=length+1; r=r-1; return true; } } } if ( (c-1) >= 0){ if (ds_cube->get_spot(s, r, (c-1)) == false) { if(path_cube->set_spot_1((s), (r), (c-1))==true){ length=length+1; c=c-1; return true; } } } if ( (c+1) < wX){ if (ds_cube->get_spot(s, r, (c+1)) == false) { if(path_cube->set_spot_1((s), (r), (c+1))==true){ length=length+1; c=c+1; return true; } } } if ( (r+1) < wY){ if (ds_cube->get_spot(s, (r+1), c) == false) { if(path_cube->set_spot_1((s), (r+1), (c))==true){ length=length+1; r=r+1; return true; } } } /*------------ SECOND PHASE OF SAME SLICE ----------------*/ if ( ((r-1) >= 0) && ((c-1) >= 0) ){ if (ds_cube->get_spot(s, (r-1), (c-1)) == false) { if(path_cube->set_spot_1((s), (r-1), (c-1))==true){ length=length+1.414; r=r-1; c=c-1; return true; } } } if ( ((r-1) >= 0) && ((c+1) < wX) ){ if (ds_cube->get_spot(s, (r-1), (c+1)) == false) { if(path_cube->set_spot_1((s), (r-1), (c+1))==true){ length=length+1.414; r=r-1; c=c+1; return true; } } } if ( ((r+1) < wY) && ((c-1) >= 0) ){ if (ds_cube->get_spot(s, (r+1), (c-1)) == false) { if(path_cube->set_spot_1((s), (r+1), (c-1))==true){ length=length+1.414; r=r+1; c=c-1; return true; } } } if ( ((r+1) < wY) && ((c+1) < wX) ){ if (ds_cube->get_spot(s, (r+1), (c+1)) == false) { if(path_cube->set_spot_1((s), (r+1), (c+1))==true){ length=length+1.414; r=r+1; c=c+1; return true; } } } if (s>0) { /*--------------- HIGHER SLICE FIRST PHASE ----------------*/ if ( (r-1) >= 0 ){ if (ds_cube->get_spot((s-1), (r-1), c) == false) { if(path_cube->set_spot_1((s-1), (r-1), (c))==true){ length=length+1.414; s=s-1; r=r-1; return true; } } } if ( (c-1) >= 0 ){ if (ds_cube->get_spot((s-1), r, (c-1)) == false) { if(path_cube->set_spot_1((s-1), (r), (c-1))==true){ length=length+1.414; s=s-1; c=c-1; return true; } } } if ( (c+1) < wX ){ if (ds_cube->get_spot((s-1), r, (c+1)) == false) { if(path_cube->set_spot_1((s-1), (r), (c+1))==true){ length=length+1.414; s=s-1; c=c+1; return true; } } } if ( (r+1) < wY ){ if (ds_cube->get_spot((s-1), (r+1), c) == false) { if(path_cube->set_spot_1((s-1), (r+1), (c))==true){ length=length+1.414; s=s-1; r=r+1; return true; } } } /*--------------- HIGHER SLICE SECOND PHASE --------------*/ if ( ((r-1) >= 0) && ((c-1) >= 0) ){ if (ds_cube->get_spot((s-1), (r-1), (c-1)) == false) { if(path_cube->set_spot_1((s-1), (r-1), (c-1))==true){; length=length+1.732; s=s-1; r=r-1; c=c-1; return true; } } } if ( ((r-1) >= 0) && ((c+1) < wX) ){ if (ds_cube->get_spot((s-1), (r-1), (c+1)) == false) { if(path_cube->set_spot_1((s-1), (r-1), (c+1))==true){ length=length+1.732; s=s-1; r=r-1; c=c+1; return true; } } } if ( ((r+1) < wY) && ((c-1) >= 0) ){ if (ds_cube->get_spot((s-1), (r+1), (c-1)) == false) { if(path_cube->set_spot_1((s-1), (r+1), (c-1))==true){ length=length+1.732; s=s-1; r=r+1; c=c-1; return true; } } } if ( ((r+1) < wY) && ((c+1) < wX) ){ if (ds_cube->get_spot((s-1), (r+1), (c+1)) == false) { if(path_cube->set_spot_1((s-1), (r+1), (c+1))==true){ length=length+1.732; s=s-1; r=r+1; c=c+1; return true; } } } /*--------- HIGHER SLICE THIRD PHASE ----------------*/ if (ds_cube->get_spot((s-1), r, c) == false) { if(path_cube->set_spot_1((s-1), (r), (c))==true){ length=length+1; s=s-1; return true; } } } return false; /*--------- DEAD PORE --------*/ }