#include "stdafx.h" #include "IsolatedVoxRemover.h" using namespace std; //constructor IsolatedVoxRemover::IsolatedVoxRemover(DS_1b_cube* input_cube) { cube = new DS_2b_cube(input_cube); wZ = input_cube->wZ; wY = input_cube->wY; wX = input_cube->wX; } //destructor IsolatedVoxRemover::~IsolatedVoxRemover() { delete cube; for (unsigned int i = 0; i < MAX_DEPTH+10; i++) delete stack[i]; delete stack; } ////recursive, if large pore can use a lot of stack space //void IsolatedVoxRemover::mark6_speed(unsigned short int zseed, unsigned short int yseed, unsigned short int xseed) //{ // if (volume >= MAX_DEPTH) // return; // // unsigned char spot_val; // // //upper spot (z-1) // int temp = zseed-1; // if ( temp >= 0 ) // { // spot_val = cube->get_spot(temp, yseed, xseed); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, temp, yseed, xseed); // mark6_speed(temp, yseed, xseed); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // //(y-1) spot // temp = yseed-1; // if ( temp >= 0) // { // spot_val = cube->get_spot(zseed, temp, xseed); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, zseed, temp, xseed); // if ((!valid_flag) && (min_Y > temp)) // min_Y = temp; // mark6_speed(zseed, temp, xseed); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // //left spot (x-1) // temp = xseed-1; // if ( temp >= 0) // { // spot_val = cube->get_spot(zseed, yseed, temp); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, zseed, yseed, temp); // if ((!valid_flag) && (min_X > temp)) // min_X = temp; // mark6_speed(zseed, yseed, temp); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // //right spot (x+1) // temp = xseed+1; // if ( temp < wX) // { // spot_val = cube->get_spot(zseed, yseed, temp); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, zseed, yseed, temp); // if ((!valid_flag) && (max_X < temp)) // max_X = temp; // mark6_speed(zseed, yseed, temp); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // //(y+1) spot // temp = yseed+1; // if ( temp < wY) // { // spot_val = cube->get_spot(zseed, temp, xseed); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, zseed, temp, xseed); // if ((!valid_flag) && (max_Y < temp)) // max_Y = temp; // mark6_speed(zseed, temp, xseed); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // //lower spot (z+1) // temp = zseed+1; // if ( temp < wZ) // { // spot_val = cube->get_spot(temp, yseed, xseed); // if (spot_val == active_val) // { // ++volume; // cube->set_spot(mark_val, temp, yseed, xseed); // if ( (!valid_flag) && (max_Z < temp) ) // max_Z = temp; // mark6_speed(temp, yseed, xseed); // } // else if ( (!valid_flag) && (spot_val == valid_val)) // { // valid_flag = true; part_flag = true; // mark_val = valid_val; // } // } // // //check if min valid volume has been reached // if ( (!valid_flag) && (volume >= min_volume) ) // { // mark_val = valid_val; // valid_flag = true; // } //} //using manual stack void IsolatedVoxRemover::mark6_stack() { unsigned short int zseed = max_Z; unsigned short int yseed = max_Y; unsigned short int xseed = max_X; stack[1]->z = max_Z; stack[1]->y = max_Y; stack[1]->x = max_X; stack[1]->next_dir = 1; unsigned char spot_val; int temp; while (stack_ptr) { if (volume >= MAX_DEPTH) { valid_flag = true; return; } zseed = stack[stack_ptr]->z; yseed = stack[stack_ptr]->y; xseed = stack[stack_ptr]->x; switch (stack[stack_ptr]->next_dir) { case 0: //done with current voxel stack_ptr--; break; case 1: //upper spot (z-1) temp = zseed - 1; if ( temp >= 0 ) { //get the location value spot_val = cube->get_spot((unsigned short)temp, yseed, xseed); if (spot_val == active_val) //if its the active value (pore or fiber) { ++volume; cube->set_spot(mark_val, (unsigned short)temp, yseed, xseed); //mark it //set the next direction to y-1 spot stack[stack_ptr]->next_dir = 2; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = (unsigned short)temp; stack[stack_ptr]->y = yseed; stack[stack_ptr]->x = xseed; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } //if found spot was previously validated from an overflow, set valid flag for current pore/fiber else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //continue on to next case } case 2: //(y-1) spot temp = yseed-1; if ( temp >= 0) { spot_val = cube->get_spot(zseed, (unsigned short)temp, xseed); if (spot_val == active_val) { ++volume; cube->set_spot(mark_val, zseed, (unsigned short)temp, xseed); if ((!valid_flag) && (min_Y > temp)) min_Y = (unsigned short)temp; //set the next direction to x-1 spot stack[stack_ptr]->next_dir = 3; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = zseed; stack[stack_ptr]->y = (unsigned short)temp; stack[stack_ptr]->x = xseed; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //continue on to next case } case 3: //left spot (x-1) temp = stack[stack_ptr]->x - 1; if ( temp >= 0) { spot_val = cube->get_spot(zseed, yseed, (unsigned short)temp); if (spot_val == active_val) { ++volume; cube->set_spot(mark_val, zseed, yseed, (unsigned short)temp); if ((!valid_flag) && (min_X > temp)) min_X = (unsigned short)temp; //set the next direction to x+1 spot stack[stack_ptr]->next_dir = 4; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = zseed; stack[stack_ptr]->y = yseed; stack[stack_ptr]->x = (unsigned short)temp; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //continue on to next case } case 4: //right spot (x+1) temp = xseed+1; if ( temp < wX) { spot_val = cube->get_spot(zseed, yseed, (unsigned short)temp); if (spot_val == active_val) { ++volume; cube->set_spot(mark_val, zseed, yseed, (unsigned short)temp); if ((!valid_flag) && (max_X < temp)) max_X = (unsigned short)temp; //set the next direction to y+1 spot stack[stack_ptr]->next_dir = 5; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = zseed; stack[stack_ptr]->y = yseed; stack[stack_ptr]->x = (unsigned short)temp; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //continue on to next case } case 5: //(y+1) spot temp = yseed+1; if ( temp < wY) { spot_val = cube->get_spot(zseed, (unsigned short)temp, xseed); if (spot_val == active_val) { ++volume; cube->set_spot(mark_val, zseed, (unsigned short)temp, xseed); if ((!valid_flag) && (max_Y < temp)) max_Y = (unsigned short)temp; //set the next direction to z+1 spot stack[stack_ptr]->next_dir = 6; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = zseed; stack[stack_ptr]->y = (unsigned short)temp; stack[stack_ptr]->x = xseed; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //continue on to next case } case 6: //lower spot (z+1) temp = zseed+1; if ( temp < wZ) { spot_val = cube->get_spot((unsigned short)temp, yseed, xseed); if (spot_val == active_val) { ++volume; cube->set_spot(mark_val, (unsigned short)temp, yseed, xseed); if ( (!valid_flag) && (max_Z < temp) ) max_Z = (unsigned short)temp; //this is last case, so next time we are done with this voxel stack[stack_ptr]->next_dir = 0; //increase stack_ptr, and fill in next spot stuff stack_ptr++; stack[stack_ptr]->z = (unsigned short)temp; stack[stack_ptr]->y = yseed; stack[stack_ptr]->x = xseed; stack[stack_ptr]->next_dir = 1; //break out, since now we'll want to process this new spot break; } else if ( (!valid_flag) && (spot_val == valid_val)) { valid_flag = true; part_flag = true; mark_val = valid_val; } //done } default: stack_ptr--; } //check if min valid volume has been reached if ( (!valid_flag) && (volume >= min_volume) ) { mark_val = valid_val; valid_flag = true; } } } void IsolatedVoxRemover::region_6(int *error_code) { //call: region_4(x,1,&err); /*Locate a black (0) region and mark it with mark_val. 4-connected.*/ register unsigned short int i, j, k; *error_code = 0; for (i = min_Z; i < wZ; i++) { for (j = ns_Y; j < wY; j++) { for (k = ns_X; k < wX; k++) { if( cube->get_spot(i, j, k) == active_val) goto pore_found; } ns_X = 0; } ns_Y = 0; } *error_code=NO_REGION; return; pore_found: max_Z = min_Z = i; ns_Y = max_Y = min_Y = j; ns_X = max_X = min_X = k; //increment volume ++volume; //mark this spot cube->set_spot(mark_val, i, j, k); mark6_stack(); //mark6_speed(i, j, k); } //call this to remove last found pore void IsolatedVoxRemover::del_reg() { register unsigned short int i, j, k; for(i = min_Z; i <= max_Z; i++) for(j = min_Y; j <= max_Y; j++) for (k = min_X; k <= max_X; k++) if(cube->get_spot(i, j, k) == mark_val) cube->set_spot(background_val, i, j, k); } //call this to validate last found pore void IsolatedVoxRemover::val_reg() { register short unsigned int i, j, k; unsigned long int cnt = 0; for(i = min_Z; i <= max_Z; i++) for(j = min_Y; j <= max_Y; j++) for (k = min_X; k <= max_X; k++) if(cube->get_spot(i, j, k) == mark_val) { cube->set_spot(valid_val, i, j, k); ++cnt; } } //returns 0 if all good, otherwise error DS_1b_cube* IsolatedVoxRemover::Run_Vox_Remover(char* out_fname, unsigned char val_to_remove, unsigned long int min_val_volume, void * dlg) { unsigned long int total_pores = 0; unsigned long int deleted_pores = 0; unsigned long int valid_pores = 0; unsigned long int part_pores = 0; int err = 0; ofstream ofile; bool to_save; if (out_fname[0] == '\0') to_save = false; else to_save = true; min_volume = min_val_volume; min_Z = ns_X = ns_Y = 0; //allocate stack MAX_DEPTH = min_volume+10000; stack = new voxel* [MAX_DEPTH+10]; if (!stack) write_to_log("isolated vox remover: could not allocate stack"); for (unsigned int i = 0; i < MAX_DEPTH+10; i++) stack[i] = new voxel; PoroMediaDialogShare* FiltDlg = (PoroMediaDialogShare*) dlg; ////progress control FiltDlg->pbar_simulation_progress.SetPos(0); FiltDlg->pbar_simulation_progress.SetStep(1); FiltDlg->pbar_simulation_progress.SetRange(0, wZ); if (val_to_remove == 0x00) //remove pores { active_val = 0; valid_val = 1; mark_val = 2; background_val = 3; } else if (val_to_remove == 0x01) //remove walls { active_val = 3; valid_val = 2; mark_val = 1; background_val = 0; } else //invalid remove val { write_to_log("Isolated Vox Remover: invalid value to remove passed (should be 0 or 1)"); return NULL; } //create the voxel.dis file if (to_save) { ofile.open(out_fname); if(!ofile.is_open()) { char log_txt[512]; sprintf(log_txt, "Isolated Vox Remover: Cannot open file %s.", out_fname); write_to_log(log_txt); return NULL; } } unsigned short int next_z = 1; /* process cube */ do { valid_flag = false; part_flag = false; stack_ptr = 1; volume = 0; region_6(&err); //locate pore if(err) break; if (!part_flag) ++total_pores; else ++part_pores; //ofile << total_pores << "\t" << volume << "\t" << min_X << "-" << max_X << "\t" << min_Y << "-" << max_Y << "\t" << min_Z << "-" << max_Z; if (val_to_remove == 0x00) mark_val = 2; else mark_val = 1; //check if volume is less than min volume, and that valid flag is NOT set if ( !valid_flag ) { del_reg(); ++deleted_pores; // ofile << "\tREMOVED\n"; } else { val_reg(); if (!part_flag) ++valid_pores; // ofile << "\n"; } if (!(total_pores % 10000)) { char log_txt[512]; sprintf(log_txt, "cur depth: %d, total: %d, removed: %d, validated: %d (part of other: %d)",min_Z,total_pores,deleted_pores,valid_pores,part_pores); write_to_log(log_txt); } //progress if (min_Z >= next_z) { FiltDlg->pbar_simulation_progress.StepIt(); ++next_z; } } while (err == 0); if (to_save) { ofile << "total found: " << total_pores << "\tremoved: " << deleted_pores << "\tvalidated: " << valid_pores << "\tpart of other: "<< part_pores << endl; ofile.close(); } DS_1b_cube * bin_cube = new DS_1b_cube(wZ, wY, wX); if (cube->to_DS_1b_cube(bin_cube, 0, wZ-1)) return NULL; return bin_cube; }