//v1.0d #include "stdafx.h" #include "bypass.h" using namespace std; void bypass1(Photoinfo* info, void * dlg_share) {//downsample directly on a FILE unsigned char n=2; //downsampling parameter (< 255) double n_cube = pow((double)n, 3); PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "by %d...", n); dlg->txt_field_progress.SetWindowText(win_txt); char* infile = info->raw_prefix; //open the slice ifstream fin(infile, ios_base::binary); if(!fin.is_open()) { char log_txt[512]; sprintf(log_txt, "DS_cube: error: cannot open input file %s", infile); write_to_log(log_txt); return; } unsigned short int Z2 = (unsigned short int)(info->sizeofz / n); //new dimensions unsigned short int Y2 = (unsigned short int)(info->sizeofy / n); //after downsampling unsigned short int X2 = (unsigned short int)(info->sizeofx / n); // //save downsampled file sprintf(win_txt, "%sdspl%d_%s(%d_%d_%d).stack.raw", dlg->field_directory_work_str, n, dlg->field_name_output_str, X2, Y2, Z2); char* outfile = win_txt; ofstream fout(outfile, ios_base::binary); if(!fout.is_open()) { char text[512]; sprintf(text,"DS_cube::write_RAW_STACK: error opening RAW stack file %s", outfile); write_to_log(text); return; } //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, Z2); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); DS_cube* cube = new DS_cube(n, info->sizeofy, info->sizeofx); //original slices cube DS_cube* cube2 = new DS_cube(1, Y2, X2); //downsampled slice for(unsigned short int z = 0; z < Z2; z++) { //read n images into memory fin.read(reinterpret_cast<char*>(&cube->memory[0]), cube->wY*cube->wX * n); { //display message in filter window sprintf(win_txt, "by %d (depth %d of %d)", n, z+1, Z2); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //downsample unsigned short int x, y; for (x=0; x < X2; x++) for (y=0; y < Y2; y++) { unsigned short int x2 = n * x; unsigned short int y2 = n * y; //convolution float total = 0; for (unsigned char i=0; i<n; i++) for (unsigned char j=0; j<n; j++) for (unsigned char k=0; k<n; k++) total += cube->data[k][y2+j][x2+i]; //save value cube2->data[0][y][x] = (unsigned char)(total/n_cube + .5); } } //save downsampled slice to disk fout.write(reinterpret_cast<char*>(cube2->memory), cube2->wZ*cube2->wY*cube2->wX); } delete cube2; delete cube; fout.close(); //downsample fin.close(); //input } void bypass2(Photoinfo* info, void * dlg_share) {//substract a FILE directly on a FILE unsigned char n=50; // # files to process at a single time (bigger is faster but takes more memory) (< 255) double n_cube; n_cube = pow((double)n, 3); PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "by %d...", n); dlg->txt_field_progress.SetWindowText(win_txt); char* infile = info->raw_prefix; //open the source slices ifstream fin(infile, ios_base::binary); if(!fin.is_open()) { char log_txt[512]; sprintf(log_txt, "DS_cube: error: cannot open input file %s", infile); write_to_log(log_txt); return; } unsigned short int Z2 = (unsigned short int)(info->sizeofz / n); //new dimensions //unsigned short int Y2 = (unsigned short int)(info->sizeofy / n); //after downsampling //unsigned short int X2 = (unsigned short int)(info->sizeofx / n); // //open file to subtract (called "[source cube name]_av") sprintf(win_txt, "%s_av", info->raw_prefix); char* infileminus = win_txt; ifstream finminus(infileminus, ios_base::binary); if(!finminus.is_open()) { char log_txt[512]; sprintf(log_txt, "DS_cube: error: cannot open input file %s", infileminus); write_to_log(log_txt); return; } //save file after minusing sprintf(win_txt, "%s_minus", info->raw_prefix); char* outfileminus = win_txt; ofstream foutminus(outfileminus, ios_base::binary); if(!foutminus.is_open()) { char text[512]; sprintf(text,"DS_cube::write_RAW_STACK: error opening RAW stack file %s", outfileminus); write_to_log(text); return; } //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, Z2); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); DS_cube* cube = new DS_cube(n, info->sizeofy, info->sizeofx); //original slices cube DS_cube* cubeminus = new DS_cube(1, cube->wY, cube->wX); //slice to subtract cube //read file to minus finminus.read(reinterpret_cast<char*>(&cubeminus->memory[0]), cubeminus->wY*cubeminus->wX*1); finminus.close(); for(unsigned short int z = 0; z < Z2; z++) { //read n images into memory fin.read(reinterpret_cast<char*>(&cube->memory[0]), cube->wY*cube->wX * n); { //display message in filter window sprintf(win_txt, "by %d (depth %d of %d)", n, z+1, Z2); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //subtracting unsigned short int x, y; for (x=0; x < cube->wX; x++) for (y=0; y < cube->wY; y++) for (unsigned char k=0; k<cube->wZ; k++) { if (cube->data[k][y][x] <= cubeminus->data[0][y][x]) cube->data[k][y][x] = 0; else cube->data[k][y][x] = cube->data[k][y][x]-cubeminus->data[0][y][x]; } } //save minused cube to disk foutminus.write(reinterpret_cast<char*>(cube->memory), cube->wZ*cube->wY*cube->wX); } delete cube; delete cubeminus; foutminus.close(); //minus fin.close(); //input } void bypass3(Photoinfo* info, void * dlg_share) {//average a FILE to a single slice ++ includes thresholding unsigned char n=50; // # files to process at a single time (bigger is faster but takes more memory) (< 255) double n_cube; n_cube = pow((double)n, 3); PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "by %d...", n); dlg->txt_field_progress.SetWindowText(win_txt); char* infile = info->raw_prefix; //open the slice ifstream fin(infile, ios_base::binary); if(!fin.is_open()) { char log_txt[512]; sprintf(log_txt, "DS_cube: error: cannot open input file %s", infile); write_to_log(log_txt); return; } unsigned short int Z2 = (unsigned short int)(info->sizeofz / n); //new dimensions //unsigned short int Y2 = (unsigned short int)(info->sizeofy / n); //after downsampling //unsigned short int X2 = (unsigned short int)(info->sizeofx / n); // //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, Z2); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); DS_cube* cube = new DS_cube(n, info->sizeofy, info->sizeofx); //original slices cube DS_ull_cube* cube_tot = new DS_ull_cube(1, cube->wY, cube->wX); //average cube //initialise av cube { for (unsigned short int j = 0; j<cube_tot->wY; j++) for (unsigned short int i = 0; i<cube_tot->wX; i++) cube_tot->data[0][j][i] = 0; } for(unsigned short int z = 0; z < Z2; z++) { //read n images into memory fin.read(reinterpret_cast<char*>(&cube->memory[0]), cube->wY*cube->wX * n); { //display message in filter window sprintf(win_txt, "by %d (depth %d of %d)", n, z+1, Z2); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //averaging unsigned short int x, y; for (x=0; x < cube->wX; x++) for (y=0; y < cube->wY; y++) for (unsigned char k=0; k<cube->wZ; k++) { if (cube->data[k][y][x] < 91) //to get rid of fibers ++ thresholding cube_tot->data[0][y][x] += cube->data[k][y][x]; else cube_tot->data[0][y][x] += 90; } } } //save results (averaging) DS_cube* cube_tot2 = new DS_cube(cube_tot->wZ, cube_tot->wY, cube_tot->wX); //average cube_tot { for (unsigned short int j = 0; j<cube_tot->wY; j++) for (unsigned short int i = 0; i<cube_tot->wX; i++) cube_tot2->data[0][j][i] = (unsigned char)(double(cube_tot->data[0][j][i])/(double)(Z2*n) + .5); } char file_out [512]; sprintf(file_out, "%s_av", dlg->field_directory_and_file_input_str); cube_tot2->write_to_RAW_STACK(file_out); delete cube_tot2; delete cube_tot; delete cube; fin.close(); //input } void bypass4(DS_cube* ds_input, void * dlg_share, char * file_name) {//histogram ds_input PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; dlg->txt_field_progress.SetWindowText("histogram..."); //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); unsigned short int height = 256; DS_cube* hist = new DS_cube(1, height+10, 260); //resulting histogram //initialise hist (image) { //white background unsigned long int k; for (k=0; k< (unsigned long int)(hist->wY*hist->wX); k++) hist->memory[k] = 255; //grayscale pattern unsigned short int i, j; for (j=height+5; j<hist->wY; j++) for (i=0; i<256; i++) hist->data[0][j][i] = (unsigned char)i; //horiz j=height; for (i=0; i<256; i++) hist->data[0][j][i] = 0; j=height+1; for (i=15; i<256; i+=16) hist->data[0][j][i] = 0; j=height+2; for (i=31; i<256; i+=32) hist->data[0][j][i] = 0; j=height+3; for (i=63; i<256; i+=64) hist->data[0][j][i] = 0; j=height+4; i=127; hist->data[0][j][i] = 0; //vert i=256; for (j=0; j<height; j++) hist->data[0][j][i] = 0; i=257; for (j=(unsigned short int)((float)height/(float)8 + .5 -1); j<height; j=j+(unsigned short int)((float)height/(float)8 + .5)) hist->data[0][j][i] = 0; i=258; for (j=(unsigned short int)((float)height/4 + .5 -1); j<height; j=j+(unsigned short int)((float)height/4 + .5)) hist->data[0][j][i] = 0; i=259; for (j=(unsigned short int)((float)height/2 + .5 -1); j<height; j=j+(unsigned short int)((float)height/2 + .5)) hist->data[0][j][i] = 0; } unsigned long long * histogram = new unsigned long long [256]; unsigned short int * histogram_display = new unsigned short int [256]; //initialise histogram (table) { unsigned short int i; for (i=0; i < 256; i++) histogram[i]=0; } //fill histogram (table) unsigned long long totnum; //of counted voxels unsigned char ignore = 1; //value to ignore [ignore] border voxels { unsigned short int i, j, k; totnum = (ds_input->wZ-2*ignore)*(ds_input->wY-2*ignore)*(ds_input->wX-2*ignore); for (k=ignore; k<ds_input->wZ-ignore; k++) { for (j=ignore; j<ds_input->wY-ignore; j++) for (i=ignore; i<ds_input->wX-ignore; i++) histogram[ ds_input->data[k][j][i] ] ++; //update progress bar dlg->pbar_simulation_progress.StepIt(); ////interrupt? //if (dlg->stop_request) // goto delvar; } } //find global max unsigned long long maxnum=0; unsigned char max_pos=0; { unsigned short int i; for (i=0; i < 256; i++) if (histogram[i] > maxnum) { maxnum = histogram[i]; max_pos = (unsigned char)i; } } //otsu thresholding: 2-levels double mu1, mu2; double sigma1, sigma2; double * histogram_double = new double [256]; { unsigned char sigma_max=0; //pi's //double * histogram_double = new double [256]; { unsigned short int i; for (i=0; i < 256; i++) histogram_double[i] = ( (double)histogram[i] / (double)maxnum ); } //sigmaB double * sigmab = new double [256]; { double w1, w2=0; unsigned short int i; for (i=0; i < 255; i++) { unsigned short int j; //1: [0,i] w1 = 0; mu1 = 0; for (j=0; j<=i; j++) { w1 += histogram_double[j]; mu1 += j * histogram_double[j]; } if (w1 == 0) goto next; mu1 /= w1; //2: [i+1, 255] w2 = 0; mu2 = 0; for (j=i+1; j<256; j++) { w2 += histogram_double[j]; mu2 += j * histogram_double[j]; } if (w2 == 0) goto next; mu2 /= w2; //sigmaB double mut = w1*mu1 + w2*mu2; sigmab[i] = w1 * pow(mu1-mut, 2) + w2 * pow(mu2-mut, 2) ; next: if ( (w1 == 0) | (w2 == 0) ) sigmab[i] = 0; } } //find gloabl max of sigmaB { double maxnum = 0; unsigned short int i; for (i=0; i < 256; i++) if (sigmab[i] > maxnum) { maxnum = sigmab[i]; sigma_max = (unsigned char)i; } } //conclude with stats { double w1, w2; unsigned short int j; //1: [0,sigma_max] w1 = 0; mu1 = 0; for (j=0; j<=sigma_max; j++) { w1 += histogram_double[j]; mu1 += j * histogram_double[j]; } mu1 /= w1; sigma1 = 0; for (j=0; j<=sigma_max; j++) sigma1 += pow(j-mu1, 2) * histogram_double[j]; sigma1 /= w1; sigma1 = pow(sigma1, .5); //2: [sigma_max+1, 255] w2 = 0; mu2 = 0; for (j=sigma_max+1; j<256; j++) { w2 += histogram_double[j]; mu2 += j * histogram_double[j]; } mu2 /= w2; sigma2 = 0; for (j=sigma_max+1; j<256; j++) sigma2 += pow(j-mu2, 2) * histogram_double[j]; sigma2 /= w2; sigma2 = pow(sigma2, .5); } } //peaks + valley finding //global max: max_pos & maxnum //local max + valley unsigned long long localmax_value=0; unsigned char localmax_pos=0; unsigned long long localmin_value; unsigned char localmin_pos; double sigma_max, sigma_loc, sigma_min; { unsigned short int i; for (i=1; i < 255; i++) //ignore 0 & 255 { //find min unsigned short int j; //unsigned char pos=i; unsigned long long value = histogram[i]; //for (j=i; j < 256; j++) if (i<max_pos) // >/< than the global max { for (j=i; j < max_pos; j++) if (value > histogram[j]) value = histogram[j]; } else for (j=max_pos; j < i; j++) if (value > histogram[j]) value = histogram[j]; //combine min&max to find local max unsigned long long val_temp; val_temp = histogram[i] - value; //find max of val_temp if (localmax_value < val_temp) { localmax_value = val_temp; localmax_pos = (unsigned char)i; } } //local min (bottom of valley between local max and gloabal max) { unsigned short int j; localmin_pos = max_pos; localmin_value = maxnum; if (localmax_pos > max_pos) { for (j=max_pos; j <= localmax_pos; j++) if (localmin_value > histogram[j]) { localmin_value = histogram[j]; localmin_pos = (unsigned char)j; } } else { for (j=localmax_pos; j <= max_pos; j++) if (localmin_value > histogram[j]) { localmin_value = histogram[j]; localmin_pos = (unsigned char)j; } } } //conclude with stats { double w; unsigned short int j; if (localmax_pos > max_pos) { //1: max_pos : [0; localmin_pos] w = 0; sigma_max = 0; for (j=0; j<=localmin_pos-1; j++) { w += histogram_double[j]; sigma_max += pow((float)j-max_pos, 2) * histogram_double[j]; } sigma_max /= w; sigma_max = pow(sigma_max, .5); //2: localmax_pos : [localmin_pos, 255] w = 0; sigma_loc = 0; for (j=localmin_pos+1; j<256; j++) { w += histogram_double[j]; sigma_loc += pow((float)j-localmax_pos, 2) * histogram_double[j]; } sigma_loc /= w; sigma_loc = pow(sigma_loc, .5); //3: localmin_pos : [max_pos localmax_pos] w = 0; sigma_min = 0; for (j=max_pos+1; j<localmax_pos; j++) { w += histogram_double[j]; sigma_min += pow((float)j-localmin_pos, 2) * histogram_double[j]; } sigma_min /= w; sigma_min = pow(sigma_min, .5); } else { //1: localmax_pos : [0; localmin_pos] w = 0; sigma_loc = 0; for (j=0; j<=localmin_pos-1; j++) { w += histogram_double[j]; sigma_loc += pow((float)j-localmax_pos, 2) * histogram_double[j]; } sigma_loc /= w; sigma_loc = pow(sigma_loc, .5); //2: max_pos : [localmin_pos, 255] w = 0; sigma_max = 0; for (j=localmin_pos+1; j<256; j++) { w += histogram_double[j]; sigma_max += pow((float)j-max_pos, 2) * histogram_double[j]; } sigma_max /= w; sigma_max = pow(sigma_max, .5); //3: localmin_pos : [localmax_pos max_pos] w = 0; sigma_min = 0; for (j=localmax_pos+1; j<max_pos; j++) { w += histogram_double[j]; sigma_min += pow((float)j-localmin_pos, 2) * histogram_double[j]; } sigma_min /= w; sigma_min = pow(sigma_min, .5); } } } //resize histogram (table) to fit display (image) double porosity128 =0; {//porosity128 -> pores = black = I<128 unsigned short int i; for (i=0; i < 128; i++) porosity128 += (double)histogram[i]; porosity128 /= totnum; } //resize histogram (table) to fit display (image) { unsigned short int i; for (i=0; i < 256; i++) histogram_display[i] = (unsigned short int)(.5 + histogram[i] * height / maxnum); } //plot hist (image) { unsigned short int i, j; //for (i=0; i < 256; i++) //full // for (j=0; j<height; j++) // if (histogram_display[i] + j >= height) // hist->data[0][j][i] = 0; //outline of histogram only / black on white //for (i=0; i < 256; i++) // for (j=0; j<height; j++) // if (histogram_display[i] + j == height) // hist->data[0][j][i] = 0; // else // hist->data[0][j][i] = 255; //outline - with dots connected / black on white i=0; hist->data[0][height - histogram_display[i]][i] = 0; for (i=1; i < 256; i++) { //numbers of dots to connect short int number = (short int)histogram_display[i] - (short int)histogram_display[i-1]; if (number == 0) hist->data[0][height - histogram_display[i]][i] = 0; else { for (j=0; j<=abs(number); j++) hist->data[0][height-histogram_display[i] + j*number/abs(number)][i -1*(int)j/abs(number)] = 0; } } //for (i=0; i < 256; i++) //outline of histogram only / white on black // for (j=0; j<height; j++) // if (histogram_display[i] + j != height) // hist->data[0][j][i] = 0; } //add grid (dots) (image) { short int i, j; unsigned char freqh = (unsigned char)(height/8); unsigned char freqv = 32; for (j=(short)(height-freqh-1); j>0; j=j-(short)(freqh)) for (i=(short)(freqv-1); i<256-freqv; i=i+(short)freqv) hist->data[0][j][i] = 192; } //save hist image char file_out [512]; //if (file_name == "") //to avoid "warning C4130: '==' : logical operation on address of string constant " if ( 0 == strcmp(file_name, "") ) sprintf(file_out, "%s_hist_%dx%d.raw", dlg->field_directory_and_file_input_str, hist->wX, hist->wY); else sprintf(file_out, "%s_hist_%dx%d.raw", file_name, hist->wX, hist->wY); hist->write_to_RAW_STACK(file_out); write_to_log(file_out); //save histogram values into file (text) { sprintf(file_out, "%s.csv", file_out); char* outfilehistogram = file_out; ofstream fout; fout.open(outfilehistogram, ios::out | ios::app); if(!fout.is_open()) { char text[512]; sprintf(text,"histogram: error opening file for saving histogram values: %s", outfilehistogram); write_to_log(text); return; } //for testing fout << (double)maxnum/(double)totnum << ", " << porosity128 << ", " \ << mu1 << ", " << sigma1 << ", " << mu2 << ", " << sigma2 << ", " \ << (unsigned short int)max_pos << ", " << sigma_max << ", " \ << (unsigned short int)localmax_pos << ", " << sigma_loc << ", " \ << (unsigned short int)localmin_pos << ", " << sigma_min << ", " \ << (double)histogram[localmin_pos]/(double)histogram[localmax_pos] << ", " \ << (double)histogram[localmin_pos]/(double)histogram[max_pos] << ", " \ << (double)histogram[localmin_pos]/(double)histogram[max_pos]*(double)histogram[localmin_pos]/(double)histogram[localmax_pos]*100 << ", " \ << 2*(double)histogram[localmin_pos]/( (double)histogram[max_pos]+(double)histogram[localmax_pos]) << ", " \ << endl << endl; //title fout << "sample, " << file_out << endl << endl; //dimensions fout << "x, " << ds_input->wX << endl; fout << "y, " << ds_input->wY << endl; fout << "z, " << ds_input->wZ << endl << endl; //ignore & total # voxels fout << "ignore, " << (unsigned short int)ignore << endl; fout << "totnum, " << totnum << endl << endl; //max fout << "max, " << maxnum << endl; fout << "max@, " << (unsigned short int)max_pos << endl; fout << "max%, " << (double)maxnum/(double)totnum << endl << endl; //porosity128 fout << "porosity128, " << porosity128 << endl << endl; //otsu results fout << "otsu 2-level" << endl; fout << "sigma_max, " << (unsigned short int)sigma_max << endl; fout << "mu1, " << mu1 << " ,sigma1, " << sigma1 << endl; fout << "mu2, " << mu2 << " ,sigma2, " << sigma2 << endl << endl; //peaks + valley results fout << "global_max@, " << (unsigned short int)max_pos << " ,value, " << histogram[max_pos] << " ,sigma_max, " << sigma_max << endl; fout << "localmax@, " << (unsigned short int)localmax_pos << " ,value, " << histogram[localmax_pos] << " ,sigma_loc, " << sigma_loc << endl; fout << "localmin@, " << (unsigned short int)localmin_pos << " ,value, " << histogram[localmin_pos] << " ,sigma_min, " << sigma_min << endl; fout << "ratio_local, " << (double)histogram[localmin_pos]/(double)histogram[localmax_pos] << endl; fout << "ratio_global, " << (double)histogram[localmin_pos]/(double)histogram[max_pos] << endl; fout << "ratio_*, " << (double)histogram[localmin_pos]/(double)histogram[max_pos]*(double)histogram[localmin_pos]/(double)histogram[localmax_pos]*100 << endl; fout << "ratio_+, " << 2*(double)histogram[localmin_pos]/( (double)histogram[max_pos]+(double)histogram[localmax_pos]) << endl << endl; //values fout << "grayscale,voxels" << endl; unsigned short int i; for (i=0; i < 256; i++) fout << i << ", " << histogram[i] << endl; fout << "****" << endl << endl; fout.close(); } // //get_threshold using minimum error // //{ // unsigned short int thresh_min_error=0; // List<Val_instance>* thresholdlist=new List<Val_instance>; // short iterationsFlag=0,imaginaryFlag=0; // register unsigned short int iterations=0; // short T, Tl = (short)lower, Tu = (short)upper; // float P1=0,P2=0; // double u1=0,u2=0,u11=0,u22=0; // double s1=0,s2=0,s11=0,s22=0; // double a=0,b=0,c1=0,d=0,D=0; // double t,t0,r1=0,r2=0,diff; // while(Tl<=Tu) // { // T=Tl; // t0=(double)T; // do // { // P1=P2=1; // u11=u22=s11=s22=0; // iterations++; // if(iterations>10000) // { // iterationsFlag=1; // break; // } // for(i=0;i<=T;i++) // P1=P1+hist[i]; // for(i=(T+1);i<=255;i++) // P2=P2+hist[i]; // for(i=0;i<=T;i++) // u11 += i*hist[i]; // for(i=(T+1); i<=255; i++) // u22 += i*hist[i]; // u1=(u11/P1); // u2=(u22/P2); // for(i=0;i<=T;i++) // s11=s11+(i-u1)*(i-u1)*hist[i]; // for(i=(T+1);i<=255;i++) // s22=s22+(i-u2)*(i-u2)*hist[i]; // s1=sqrt(s11/P1); // s2=sqrt(s22/P2); // if ((s1 ==0) || (s2 ==0)) // break; // a=(1/(s1*s1)-1/(s2*s2)); // b=-2*(u1/(s1*s1)-u2/(s2*s2)); // c1=(u1*u1)/(s1*s1)-(u2*u2)/(s2*s2); // c1=c1+2*(log(s1)-log(s2))-2*(log(P1)-log(P2)); // d=(b*b-4*a*c1); // if(d < 0) // { // imaginaryFlag = 1; // break; // } // D=sqrt(d); // r1=(-b+D)/(2*a); // r2=(-b-D)/(2*a); // if(r1>u1 && r1<u2) // t=r1; // else // t=r2; // if(t0>t) // diff=t0-t; // else // diff=t-t0; // t0 = t; // T = (short)t; // //if ((t-(double)T) > 0.5) //round to closest int // // T ++; //no, because [>] threshold in binarise.cpp // } while (diff >= 0.01); // if(iterationsFlag==1) // { // char log_txt[512]; // sprintf(log_txt, "Threshold: iterations > 10000 for Tl=%d", Tl); // write_to_log(log_txt); // iterationsFlag=0; // } // else if(imaginaryFlag==1) // { // char log_txt[512]; // sprintf(log_txt, "Threshold: b2-4ac<0 for Tl=%d", Tl); // write_to_log(log_txt); // imaginaryFlag=0; // } // else //fill list // { // if ((T>=0)&&(T<=255)) // { // bool is_add=false; // for(i=1; i<=thresholdlist->getNumber(); i++) // if(thresholdlist->getWhich(i)->value==T) // { // thresholdlist->getWhich(i)->instance++; // is_add=true; // break; // } // if(!is_add) // { // Val_instance* newthreshold=new Val_instance; // newthreshold->value=T; // newthreshold->instance=1; // thresholdlist->Add(*newthreshold); // } // } // } // Tl++; // iterations=0; // } // if (thresholdlist->getNumber() > 0) // { // result.is_get=true; // int mostcommon=0; // for(i=1; i<=thresholdlist->getNumber() ;i++) //ignore value >upper and <lower // if ((thresholdlist->getWhich(i)->instance > mostcommon)&&(thresholdlist->getWhich(i)->value <=upper)&&(thresholdlist->getWhich(i)->value >=lower)) // { // mostcommon = thresholdlist->getWhich(i)->instance; // thresh_min_error = (unsigned short)thresholdlist->getWhich(i)->value; // } // } // delete thresholdlist; // { // char log_txt[512]; // sprintf(log_txt, "Threshold: based on minimum error thresholding: %d", thresh_min_error); // write_to_log(log_txt); // } //} dlg->txt_field_progress.SetWindowText("histogram...done!"); //delvar: delete [] histogram_double; delete [] histogram; delete [] histogram_display; delete hist; } void bypass5(DS_cube* ds_input, void * dlg_share) {//unwrap an ORIGINAL ENTIRE RAW image PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; dlg->txt_field_progress.SetWindowText("unwrap..."); //find center of image unsigned short int centerx, centery; { unsigned short int i=0, j, k; k = 0; //-> use first slice only unsigned short int westx, westy1, westy2, eastx, easty1, easty2; unsigned short int northy, northx1, northx2, southy, southx1, southx2; unsigned char background = ds_input->data[0][0][0]; //define background color as top left corner voxel value for (j=0; j<ds_input->wY; j++) for (i=0; i<ds_input->wX; i++) if ( ds_input->data[k][j][i] != background ) goto out_north; out_north: northy = j; northx1 = i; for (i=ds_input->wX-1; i>0; i--) if ( ds_input->data[k][j][i] != background ) break; northx2 = i; for (j=ds_input->wY-1; j>0; j--) for (i=0; i<ds_input->wX; i++) if ( ds_input->data[k][j][i] != background ) goto out_south; out_south: southy = j; southx1 = i; for (i=ds_input->wX-1; i>0; i--) if ( ds_input->data[k][j][i] != background ) break; southx2 = i; for (i=0; i<ds_input->wX; i++) for (j=0; j<ds_input->wY; j++) if ( ds_input->data[k][j][i] != background ) goto out_west; out_west: westx = i; westy1 = j; for (j=ds_input->wY-1; j>0; j--) if ( ds_input->data[k][j][i] != background ) break; westy2 = j; for (i=ds_input->wX-1; i>0; i--) for (j=0; j<ds_input->wY; j++) if ( ds_input->data[k][j][i] != background ) goto out_east; out_east: eastx = i; easty1 = j; for (j=ds_input->wY-1; j>0; j--) if ( ds_input->data[k][j][i] != background ) break; easty2 = j; centerx = (northx2 - northx1)/2 + northx1; centerx = (southx2 - southx1)/2 + southx1; //should be identical centery = (westy2 - westy1)/2 + westy1; centery = (easty2 - easty1)/2 + easty1; //should be identical } //find width & height unsigned short int width = centerx; if (width > centery) width = centery; if ( width > (ds_input->wY - centery) ) width = ds_input->wY - centery; if ( width > (ds_input->wX - centerx) ) width = ds_input->wX - centerx; unsigned short int height; { const double pi = 3.1415926535897932384626433832795; height = (unsigned short int)(2*pi*width +.5); } DS_cube* unwrap = new DS_cube(ds_input->wZ, height, width); //unwrapped stack //unwrap { unsigned short int i, j, k; const double pi = 3.1415926535897932384626433832795; for (k=0; k < ds_input->wZ; k++) for (j=0; j < height; j++) for (i=0; i<width; i++) { double theta = 2*pi * double(j)/double(height); unwrap->data[k][j][i] = ds_input->data[k][centery+(short int)(.5+i*cos(theta))][centerx+(short int)(.5+i*sin(theta))] ; } } //save char file_out [512]; sprintf(file_out, "%s_unwrap(%d_%d_%d).stack.raw", dlg->field_directory_and_file_input_str, unwrap->wX, unwrap->wY, unwrap->wZ); unwrap->write_to_RAW_STACK(file_out); delete unwrap; } //void bypass(Photoinfo* info, void * dlg_share) //{ // unsigned char n=50; //downsampling parameter // double n_cube = pow((double)n, 3); // // PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; // char win_txt [512]; // sprintf(win_txt, "by %d...", n); // dlg->txt_field_progress.SetWindowText(win_txt); // // char* infile = info->raw_prefix; // //open the slice // ifstream fin(infile, ios_base::binary); // if(!fin.is_open()) // { // char log_txt[512]; // sprintf(log_txt, "DS_cube: error: cannot open input file %s", infile); // write_to_log(log_txt); // return; // } // // unsigned short int Z2 = (unsigned short int)(info->sizeofz / n); //new dimensions // unsigned short int Y2 = (unsigned short int)(info->sizeofy / n); //after downsampling // unsigned short int X2 = (unsigned short int)(info->sizeofx / n); // // // //open file to subtract // sprintf(win_txt, "%s_av", info->raw_prefix); // char* infileminus = win_txt; // ifstream finminus(infileminus, ios_base::binary); // if(!finminus.is_open()) // { // char log_txt[512]; // sprintf(log_txt, "DS_cube: error: cannot open input file %s", infileminus); // write_to_log(log_txt); // return; // } // // ////save downsampled file // //sprintf(win_txt, "%sdspl%d_%s(%d_%d_%d).stack.raw", dlg->field_directory_work_str, n, dlg->field_name_output_str, X2, Y2, Z2); // //char* outfile = win_txt; // //ofstream fout(outfile, ios_base::binary); // //if(!fout.is_open()) // //{ // // char text[512]; // // sprintf(text,"DS_cube::write_RAW_STACK: error opening RAW stack file %s", outfile); // // write_to_log(text); // // return; // //} // // //save file after minusing // sprintf(win_txt, "%s_minus", info->raw_prefix); // char* outfileminus = win_txt; // ofstream foutminus(outfileminus, ios_base::binary); // if(!foutminus.is_open()) // { // char text[512]; // sprintf(text,"DS_cube::write_RAW_STACK: error opening RAW stack file %s", outfileminus); // write_to_log(text); // return; // } // // //setup the progress bar parameters // dlg->pbar_simulation_progress.SetRange(0, Z2); // dlg->pbar_simulation_progress.SetPos(0); // dlg->pbar_simulation_progress.SetStep(1); // // DS_cube* cube = new DS_cube(n, info->sizeofy, info->sizeofx); //original slices cube // DS_cube* cubeminus = new DS_cube(1, cube->wY, cube->wX); //slice to subtract cube // DS_cube* cube2 = new DS_cube(1, Y2, X2); //downsampled slice // DS_ull_cube* cube_tot = new DS_ull_cube(1, cube->wY, cube->wX); //average cube // // //read file to minus // finminus.read(reinterpret_cast<char*>(&cubeminus->memory[0]), cubeminus->wY*cubeminus->wX*1); // finminus.close(); // // ////initialise av cube // //{ // // for (unsigned short int j = 0; j<cube_tot->wY; j++) // // for (unsigned short int i = 0; i<cube_tot->wX; i++) // // cube_tot->data[0][j][i] = 0; // //} // // for(unsigned short int z = 0; z < Z2; z++) // { // //read n images into memory // //fin.clear(); // //fin.seekg(z*cube->wY*cube->wX*n); // fin.read(reinterpret_cast<char*>(&cube->memory[0]), cube->wY*cube->wX * n); // // { // //display message in filter window // sprintf(win_txt, "by %d (depth %d of %d)", n, z+1, Z2); // dlg->txt_field_progress.SetWindowText(win_txt); // //update progress bar // dlg->pbar_simulation_progress.StepIt(); // // //subtracting // unsigned short int x, y; // for (x=0; x < cube->wX; x++) // for (y=0; y < cube->wY; y++) // for (unsigned char k=0; k<cube->wZ; k++) // { // if (cube->data[k][y][x] <= cubeminus->data[0][y][x]) // cube->data[k][y][x] = 0; // else // cube->data[k][y][x] -= cubeminus->data[0][y][x]; // } // // ////averaging // // unsigned short int x, y; // // for (x=0; x < cube->wX; x++) // // for (y=0; y < cube->wY; y++) // // for (unsigned char k=0; k<cube->wZ; k++) // // { // // if (cube->data[k][y][x] < 91) //to get rid of fibers // // cube_tot->data[0][y][x] += cube->data[k][y][x]; // // else // // cube_tot->data[0][y][x] += 90; // // } // ////downsample // //unsigned short int x, y; // //for (x=0; x < X2; x++) // // for (y=0; y < Y2; y++) // // { // // unsigned short int x2 = n * x; // // unsigned short int y2 = n * y; // // //convolution // // float total = 0; // // for (unsigned char i=0; i<n; i++) // // for (unsigned char j=0; j<n; j++) // // for (k=0; k<n; k++) // // total += cube->data[k][y2+j][x2+i]; // // //save value // // cube2->data[0][y][x] = (unsigned char)(total/n_cube + .5); // // } // } // //save downsampled slice to disk // //fout.write(reinterpret_cast<char*>(cube2->memory), cube2->wZ*cube2->wY*cube2->wX); // //save minused cube to disk // foutminus.write(reinterpret_cast<char*>(cube->memory), cube->wZ*cube->wY*cube->wX); // } // // ////save results (averaging) // //DS_cube* cube_tot2 = new DS_cube(cube_tot->wZ, cube_tot->wY, cube_tot->wX); // ////average cube_tot // //{ // // for (unsigned short int j = 0; j<cube_tot->wY; j++) // // for (unsigned short int i = 0; i<cube_tot->wX; i++) // // cube_tot2->data[0][j][i] = (unsigned char)(double(cube_tot->data[0][j][i])/(double)(Z2*n) + .5); // //} // // //char file_out [512]; // //sprintf(file_out, "%s_av", dlg->field_directory_and_file_input_str); // //cube_tot2->write_to_RAW_STACK(file_out); // //delete cube_tot2; // //delete cube_tot; // // delete cube2; // delete cube; // delete cubeminus; // // //fout.close(); //downsample // foutminus.close(); //minus // fin.close(); //input //}