//v1.2a #include "convolve_filters.h" #include "StdAfx.h" #include <algorithm> using namespace std; unsigned char median_conv(unsigned char input [], unsigned char l) { //l size of input (# elements) sort(input, input + l); //built-in sort algorithm return input[l/2]; //pure median } unsigned char median_w_averaging_conv(unsigned char input [], unsigned char l) //, unsigned char n_med) { //n # points around the median used for averaging //l size of input (# elements) sort(input, input + l); //built-in sort algorithm register unsigned short int sum = 0; unsigned char n_med=0; n_med = (unsigned char)(l/3); //disregard 1/3rd of array on each side for (register unsigned short int i=n_med; i<l-n_med; i++) sum = sum + (unsigned short)input[i]; return (unsigned char)( (double)sum/(double)(l-2*n_med) + 0.5 ); //+0.5 -> round to closest int } bool convolve_filter(DS_cube* ds_input, void * dlg_share) { //perform Convolution Filter //display message in filter window PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "Convolution Filter ..."); dlg->txt_field_progress.SetWindowText(win_txt); //determine mask size unsigned short int mask_size=0; float radius = dlg->filter_convolve_radius; mask_size = 2 * (short int)(radius) + 1 ; if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) return false; //sample too small //allocate memory for the mask filter DS_f_cube* mask_filter = new DS_f_cube(mask_size, mask_size, mask_size); //set values to zero unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) mask_filter->data[k][j][i] = 0.0; //fill mask filter float radius_sq = radius*radius; unsigned short int mask_radius=0; mask_radius = (short int)(radius); //type of filling switch (dlg->filter_convolve_type) { case 2: { //uniform (MEAN) for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if ( pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2) < radius_sq) mask_filter->data[k][j][i] = 1; break; } case 3: { //distance: LINEAR for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) { register float distance_sq = pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2); if ( distance_sq < radius_sq) mask_filter->data[k][j][i] = 1 - sqrt(distance_sq / radius_sq) / 2; //linear } break; } case 4: { //distance: EXP float sigma = dlg->filter_convolve_radius; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) { register float distance_sq = pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2); if ( distance_sq < radius_sq) mask_filter->data[k][j][i] = exp((float)(-.5)*(distance_sq/sigma/sigma)); //exp } break; } case 5: { //distance: SQ for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) { register float distance_sq = pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2); if ( distance_sq < radius_sq) mask_filter->data[k][j][i] = 1/pow(distance_sq, 1/6); //square } break; } } //normalize filter { register float total = 0.0; //add up for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) total += mask_filter->data[k][j][i]; //normalize for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) mask_filter->data[k][j][i] /= total; } //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ - mask_size + 1); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //allocate mem for temporary slice buffer DS_cube* slice_buffer = new DS_cube(mask_radius + 1, ds_input->wY, ds_input->wX); //calculate buffer results for the first (mask_size - 1) slices for (unsigned short int depth = 0; depth <= mask_radius; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter ... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[depth][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); unsigned short int x,y; for (x=0; x <= ds_input->wX - mask_size; x++) for (y=0; y <= ds_input->wY - mask_size; y++) { //convolution register float total = 0.0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) total += mask_filter->data[k][j][i] * ds_input->data[depth+k][y+j][x+i]; //save value slice_buffer->data[depth][y+mask_radius][x+mask_radius] = (unsigned char)(total+0.5); } } //filter for remaining slices for (unsigned short int depth = mask_radius+1; depth <= ds_input->wZ - mask_size; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter ... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put 1st line of slice_buffer into ds_input @ depth-1 memcpy(&ds_input->data[depth - 1][0][0], &slice_buffer->data[0][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //rotate buffer (copy everything from the second slice on the first slice) //unsigned char **sbtemp; //sbtemp = slice_buffer->data[0]; //for (k=0; k < mask_radius ; k++) // slice_buffer->data[k] = slice_buffer->data[k+1]; //slice_buffer->data[mask_radius] = sbtemp; memcpy(&slice_buffer->data[0][0][0], &slice_buffer->data[1][0][0], mask_radius * ds_input->wY * ds_input->wX * sizeof(unsigned char)); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[mask_radius][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //save new results to last line of buffer unsigned short int x,y; for (x=0; x<=ds_input->wX - mask_size; x++) for (y=0; y<=ds_input->wY - mask_size; y++) { //convolution register float total = 0.0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) total += mask_filter->data[k][j][i] * ds_input->data[depth+k][y+j][x+i]; //save value slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = (unsigned char)(total+0.5); } //if last slice, save results into ds if ( depth == ds_input->wZ - mask_size) { //put whole slice_buffer into ds_input (remaining last slices) memcpy(&ds_input->data[depth][0][0], &slice_buffer->data[0][0][0], (mask_radius + 1) * ds_input->wY * ds_input->wX * sizeof(unsigned char)); } if (dlg->stop_request) { delete slice_buffer; delete mask_filter; return false; } } //clear variables //delvar: delete slice_buffer; delete mask_filter; return true; } bool convolve_filter_median(DS_cube* ds_input, void * dlg_share) { //perform Convolution Filter with MEDIAN //display message in filter window PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "Convolution Filter ..."); dlg->txt_field_progress.SetWindowText(win_txt); //pointer to median function unsigned char (*pt2median) (unsigned char [], unsigned char); switch (dlg->filter_convolve_type) { case 0: pt2median = median_conv; break; default: //to avoid compiler warning message **could be a problem if the drop down box has more values** //case 1: pt2median = median_w_averaging_conv; break; } //determine mask size unsigned short int mask_size=0; float radius = dlg->filter_convolve_radius; mask_size = 2 * (short int)(radius) + 1 ; if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) return false; //sample too small //allocate memory for the mask filter DS_usi_cube* mask_filter = new DS_usi_cube(mask_size, mask_size, mask_size); //set values to zero unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) mask_filter->data[k][j][i] = 0; //fill mask filter float radius_sq = radius*radius; unsigned short int mask_radius = 0; mask_radius = (short int)(radius); unsigned char median_count = 0; //biggest cube ~ 6*6*6 //uniform filling (for median) for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if ( pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2) < radius_sq) { mask_filter->data[k][j][i] = 1; median_count ++; } //no normalization //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ - mask_size + 1); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //allocate mem for temporary slice buffer DS_cube* slice_buffer = new DS_cube(mask_radius + 1, ds_input->wY, ds_input->wX); //mem for temporary buffer for median determination unsigned char* values = new unsigned char [median_count]; //calculate buffer results for the first (mask_size - 1) slices for (unsigned short int depth = 0; depth <= mask_radius; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter ... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[depth][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); unsigned short int x,y; for (x=0; x<=ds_input->wX - mask_size; x++) for (y=0; y<=ds_input->wY - mask_size; y++) { //convolution register float total; total = 0; register unsigned short int ind_val = 0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if (mask_filter->data[k][j][i]) { values[ind_val] = ds_input->data[depth+k][y+j][x+i]; ind_val ++; } //save value slice_buffer->data[depth][y+mask_radius][x+mask_radius] = pt2median(values, median_count); //median(values, median_count); } } //filter for remaining slices for (unsigned short int depth = mask_radius+1; depth <= ds_input->wZ - mask_size; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter ... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put 1st line of slice_buffer into ds_input @ depth-mask_radius memcpy(&ds_input->data[depth - 1][0][0], &slice_buffer->data[0][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //rotate buffer memcpy(&slice_buffer->data[0][0][0], &slice_buffer->data[1][0][0], mask_radius * ds_input->wY * ds_input->wX * sizeof(unsigned char)); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[mask_radius][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //save new results to last line of buffer unsigned short int x,y; for (x=0; x<=ds_input->wX - mask_size; x++) for (y=0; y<=ds_input->wY - mask_size; y++) { //convolution register float total; total = 0; register unsigned short int ind_val = 0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if (mask_filter->data[k][j][i]) { values[ind_val] = ds_input->data[depth+k][y+j][x+i]; ind_val ++; } //save value slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = pt2median(values, median_count); //median(values, median_count); } //if last slice, save results into ds if ( depth == ds_input->wZ - mask_size) { //put whole slice_buffer into ds_input (remaining last slices) memcpy(&ds_input->data[depth][0][0], &slice_buffer->data[0][0][0], (mask_radius + 1) * ds_input->wY * ds_input->wX * sizeof(unsigned char)); } if (dlg->stop_request) { delete slice_buffer; delete mask_filter; delete [] values; return false; } } //clear variables //delvar: delete slice_buffer; delete mask_filter; delete [] values; return true; } bool convolve_threshold(DS_cube* ds_input, void * dlg_share) { //perform Convolution Filter for thresholding //display message in filter window PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "Convolution Filter - Threshold..."); dlg->txt_field_progress.SetWindowText(win_txt); //init parameters unsigned char min = 107; unsigned char max = 148; //determine mask size float radius = dlg->filter_convolve_radius; unsigned short int mask_radius = 0; mask_radius = (short int)(radius); unsigned short int mask_size=0; mask_size = 2 * mask_radius + 1 ; if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) return false; //sample too small //allocate memory for the mask filter DS_usi_cube* mask_filter = new DS_usi_cube(mask_size, mask_size, mask_size); //set values to zero unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) mask_filter->data[k][j][i] = 0; //fill mask filter float radius_sq = radius*radius; mask_radius = (short int)(radius); //uniform filling (for median) for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if ( pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2) < radius_sq) { mask_filter->data[k][j][i] = 1; } //no normalization //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ - mask_size + 1); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //allocate mem for temporary slice buffer DS_cube* slice_buffer = new DS_cube(mask_radius + 1, ds_input->wY, ds_input->wX); //calculate buffer results for the first (mask_size - 1) slices for (unsigned short int depth = 0; depth <= mask_radius; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter - Threshold... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[depth][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); register unsigned short int x,y; for (x=0; x<=ds_input->wX - mask_size; x++) for (y=0; y<=ds_input->wY - mask_size; y++) if (ds_input->data[depth][y+mask_radius][x+mask_radius] > max) slice_buffer->data[depth][y+mask_radius][x+mask_radius] = 255; else if (ds_input->data[depth][y+mask_radius][x+mask_radius] < min) slice_buffer->data[depth][y+mask_radius][x+mask_radius] = 0; else { //convolution register unsigned char above = 0; register unsigned char below = 0; register unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if (mask_filter->data[k][j][i]) { if (ds_input->data[depth+k][y+j][x+i] > max) above ++; else if (ds_input->data[depth+k][y+j][x+i] < min) below ++; } //save value if (below > above) slice_buffer->data[depth][y+mask_radius][x+mask_radius] = 0; else if (below < above) slice_buffer->data[depth][y+mask_radius][x+mask_radius] = 255; } } //filter for remaining slices for (unsigned short int depth = mask_radius+1; depth <= ds_input->wZ - mask_size; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter ... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put 1st line of slice_buffer into ds_input @ depth-mask_radius memcpy(&ds_input->data[depth - 1][0][0], &slice_buffer->data[0][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //rotate buffer memcpy(&slice_buffer->data[0][0][0], &slice_buffer->data[1][0][0], mask_radius * ds_input->wY * ds_input->wX * sizeof(unsigned char)); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[mask_radius][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //save new results to last line of buffer register unsigned short int x,y; for (x=0; x<=ds_input->wX - mask_size; x++) for (y=0; y<=ds_input->wY - mask_size; y++) if (ds_input->data[depth][y+mask_radius][x+mask_radius] > max) slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = 255; else if (ds_input->data[depth][y+mask_radius][x+mask_radius] < min) slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = 0; else { //convolution register unsigned char above = 0; register unsigned char below = 0; register unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if (mask_filter->data[k][j][i]) { if (ds_input->data[depth+k][y+j][x+i] > max) above ++; else if (ds_input->data[depth+k][y+j][x+i] < min) below ++; } //save value if (below > above) slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = 0; else if (below < above) slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = 255; } //if last slice, save results into ds if ( depth == ds_input->wZ - mask_size ) { //put whole slice_buffer into ds_input (remaining last slices) memcpy(&ds_input->data[depth][0][0], &slice_buffer->data[0][0][0], (mask_radius + 1) * ds_input->wY * ds_input->wX * sizeof(unsigned char)); } if (dlg->stop_request) { delete slice_buffer; delete mask_filter; return false; } } //clear variables //delvar: delete slice_buffer; delete mask_filter; return true; } float difcoef1(float in, float K) { return ( in/(1+pow(in/K, 2)) ); //noise removal } float difcoefe2(float in, float K) { return ( in*exp( -pow (in/K, 2) ) ); //high contrast preservation } float difcoefe4(float in, float K) { return ( in*exp( -pow (in/K, 4) ) ); //high contrast preservation } float difcoefe20(float in, float K) { return ( in*exp( -pow (in/K, 20) ) ); //high contrast preservation } bool convolve_anisotropic_diffusion(DS_cube* ds_input, void * dlg_share) {//perform Convolution Filter for anisotropic_diffusion //display message in filter window PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "Convolution Filter - anisotropic_diffusion..."); dlg->txt_field_progress.SetWindowText(win_txt); //goto delvar; //for debug purpose //histogram original cube bypass4(ds_input, dlg, ""); //histogram only <-> (user selects K=1 in GUI) if (dlg->filter_aniso_K == 1) goto delvar; //pointer to pt2diff_coeff float(*pt2difcoef) (float, float); //aniso filter list box 0 = 1/[1+]^2, 1 = exp^2, 2 = exp^4, 3 = exp^n switch (dlg->filter_aniso_diff_coef) { case 0: pt2difcoef = difcoef1; break; case 1: pt2difcoef = difcoefe2; break; case 2: pt2difcoef = difcoefe4; break; default: //case 3: pt2difcoef = difcoefe20; break; } //counting variables for convergence determination long double total_change; long double total_change_abs; long double total_change_sq; long double old_total_change_sq=0; unsigned long long threshold_change; long long totreal_change; unsigned long long totreal_change_abs, totreal_change_sq; //decimal ; //float K=14; //for exp^4 -> 14 //float tau=1; float K=dlg->filter_aniso_K; float tau=dlg->filter_aniso_TAU; //look-up table float * difcoef_lut = new float [512]; //fill table for non variable K if (!dlg->filter_aniso_K_var) { short int i; for (i=-255; i<256; i++) difcoef_lut[255+i] = pt2difcoef(i, K); } //determine mask size float radius = dlg->filter_aniso_radius; unsigned short int mask_radius = 0; mask_radius = (short int)(radius); unsigned short int mask_size=0; mask_size = 2 * mask_radius + 1 ; if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) return false; //sample too small //allocate memory for the mask filter DS_f_cube* mask_filter = new DS_f_cube(mask_size, mask_size, mask_size); //set values to zero unsigned char i,j,k; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) mask_filter->data[k][j][i] = 0.0; //fill mask filter float radius_sq = radius*radius; mask_radius = (short int)(radius); //uniform filling (for median) for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if ( pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2) < radius_sq ) { mask_filter->data[k][j][i] = 1; } //normalize filter float totalfilter = 0; //add up -> totalfilter for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) totalfilter += mask_filter->data[k][j][i]; ////normalize //for (i=0; i<mask_size; i++) // for (j=0; j<mask_size; j++) // for (k=0; k<mask_size; k++) // mask_filter->data[k][j][i] /= total; //save original cube middle SLICE to file { //file name char file_out_slice [512]; char int_file_name [512]; sprintf(int_file_name, "%s//%s", dlg->field_directory_work_str, dlg->field_name_output_str); sprintf(file_out_slice, "%s_K%d_i#0(%d_%d_1).stack.raw", int_file_name, \ (unsigned char)dlg->filter_aniso_K, ds_input->wX, ds_input->wY); dlg->txt_field_progress.SetWindowText("Saving anisotropic diffusion filter result..."); //temporary slice DS_slice* dslice = new DS_slice(ds_input->wY, ds_input->wX); //determine depth k int k = 0; k = (int)(ds_input->wZ / 2); //extract temporary slice @ depth k ds_input->to_DS_slice(dslice, (unsigned short)k); dslice->write_to_RAW(file_out_slice); delete dslice; } //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ - mask_size + 1); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //allocate mem for temporary slice buffer DS_cube* slice_buffer = new DS_cube(mask_radius + 1, ds_input->wY, ds_input->wX); for (unsigned char iters = 1; iters <= dlg->filter_aniso_iterations; iters++) { //initialise counting variables for convergence determination total_change = 0; total_change_sq = 0; total_change_abs = 0; threshold_change = 0; totreal_change = 0; totreal_change_abs = 0; totreal_change_sq = 0; //for variable K : if (dlg->filter_aniso_K_var && ( iters < 4 ) ) { //K = (int)( 14 + dlg->filter_aniso_K*exp(-(float)(iters)/2)); switch (iters) { case 1: K = 60; break; case 2: K = 30; break; case 3: K = dlg->filter_aniso_K; break; default:; } //update LUT: fill table //if ( iters < 4 ) { short int i; for (i=-255; i<256; i++) difcoef_lut[255+i] = pt2difcoef(i, K); } } //calculate buffer results for the first (mask_size - 1) slices for (unsigned short int depth = 0; depth <= mask_radius; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter - aniso... iter#%d (depth %d of %d)", iters, depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[depth][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); unsigned short int x,y; for (x=0; x <= ds_input->wX - mask_size; x++) for (y=0; y <= ds_input->wY - mask_size; y++) { //convolution - bulk float total = 0; //unsigned short int b = 0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) //if (mask_filter->data[k][j][i]) { //unsigned short int a; //total += mask_filter->data[k][j][i] * pt2difcoef((ds_input->data[depth+k][y+j][x+i] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]), K); total += mask_filter->data[k][j][i] * difcoef_lut[ 255+ds_input->data[depth+k][y+j][x+i] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] ]; //b += a; } //save value total *= tau/totalfilter; total_change += total; total_change_abs += abs(total); total_change_sq += total*total; slice_buffer->data[depth][y+mask_radius][x+mask_radius] = (unsigned char)(ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] + total + .5); //round threshold_change += ( (slice_buffer->data[depth][y+mask_radius][x+mask_radius] >> 7) \ ^ (ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] >> 7 ) ); totreal_change += slice_buffer->data[depth][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]; totreal_change_abs += abs(slice_buffer->data[depth][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]); totreal_change_sq += (unsigned long long)pow((float)(slice_buffer->data[depth][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]), 2); } //convolution - sides/main borders (leave corners alone) //x=0; } //filter for remaining slices for (unsigned short int depth = mask_radius+1; depth <= ds_input->wZ - mask_size; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter - aniso... iter#%d (depth %d of %d)", iters, depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //put 1st line of slice_buffer into ds_input @ depth-mask_radius memcpy(&ds_input->data[depth - 1][0][0], &slice_buffer->data[0][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //rotate buffer memcpy(&slice_buffer->data[0][0][0], &slice_buffer->data[1][0][0], mask_radius * ds_input->wY * ds_input->wX * sizeof(unsigned char)); //put ds_input @ slice = depth+mask_radius into slice_buffer (to copy unprocessed borders too) memcpy(&slice_buffer->data[mask_radius][0][0], &ds_input->data[depth + mask_radius][0][0], ds_input->wY * ds_input->wX * sizeof(unsigned char)); //save new results to last line of buffer unsigned short int x,y; for (x=0; x <= ds_input->wX - mask_size; x++) for (y=0; y <= ds_input->wY - mask_size; y++) { //convolution - bulk float total = 0.0; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) total += mask_filter->data[k][j][i] * difcoef_lut[ 255+ds_input->data[depth+k][y+j][x+i] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] ]; //total += mask_filter->data[k][j][i] * pt2difcoef(ds_input->data[depth+k][y+j][x+i] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius], K); //save value total *= tau/totalfilter; total_change += total; total_change_abs += abs(total); total_change_sq += total*total; slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] = (unsigned char)(ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] + total +.5); threshold_change += ( (ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius] >> 7) \ ^ (slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] >> 7) ); totreal_change += slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]; totreal_change_abs += abs(slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]); totreal_change_sq += (unsigned long long)pow((float)(slice_buffer->data[mask_radius][y+mask_radius][x+mask_radius] - ds_input->data[depth+mask_radius][y+mask_radius][x+mask_radius]), 2); } //if last slice, save results into ds if ( depth == ds_input->wZ - mask_size ) { //put whole slice_buffer into ds_input (remaining last slices) memcpy(&ds_input->data[depth][0][0], &slice_buffer->data[0][0][0], (mask_radius + 1) * ds_input->wY * ds_input->wX * sizeof(unsigned char)); } if (dlg->stop_request) goto delvar; } //save results to file + histogram { char file_out [512]; sprintf(file_out, "%s//%s(%d_%d_%d).stack.raw", dlg->field_directory_work_str, dlg->field_name_output_str, ds_input->wX, ds_input->wY, ds_input->wZ); sprintf(file_out, "%s_K%d_i#%d", file_out, (unsigned short int)(K+.5), iters); //save filtered cube to file //if ( (((iters % 2 == 0) && (iters < 40)) || (iters % 5 == 0) || (iters < 6) ) && dlg->save_intermediate_results) //{ // dlg->txt_field_progress.SetWindowText("Saving anisotropic diffusion filter result..."); // ds_input->write_to_RAW_STACK(file_out); //} //save filtered cube middle SLICE to file (little hard drive usage) { //file name char file_out_slice [512]; char int_file_name [512]; sprintf(int_file_name, "%s//%s", dlg->field_directory_work_str, dlg->field_name_output_str); sprintf(file_out_slice, "%s_K%d_i#%d(%d_%d_1).stack.raw", int_file_name, \ (unsigned char)dlg->filter_aniso_K, iters, ds_input->wX, ds_input->wY); dlg->txt_field_progress.SetWindowText("Saving anisotropic diffusion filter result..."); //temporary slice DS_slice* dslice = new DS_slice(ds_input->wY, ds_input->wX); //determine depth k int k = 0; k = (int)(ds_input->wZ / 2); //extract temporary slice @ depth k ds_input->to_DS_slice(dslice, (unsigned short)k); dslice->write_to_RAW(file_out_slice); delete dslice; } //histogram result bypass4(ds_input, dlg, file_out); //write sample name to log //write_to_log(file_out); sprintf(file_out, "total_change, %Lf, total_change_abs, %Lf, total_change_sq, %Lf, threshold_change, %I64u, \ totreal_change, %I64d, totreal_change_abs, %I64u, totreal_change_sq, %I64u", \ total_change, total_change_abs, total_change_sq, threshold_change, \ totreal_change, totreal_change_abs, totreal_change_sq); write_to_log(file_out); } //auto-stop @ 5% (delta=(old-new)/new)@total_change_sq if ( (iters > 4) && ( old_total_change_sq/total_change_sq < 1.05) ) { dlg->filter_aniso_iterations = iters; write_to_log("# iterations = %d", iters); goto delvar; } old_total_change_sq = total_change_sq; //contrast boost if (dlg->filter_aniso_conboost) { unsigned char below, above, delta; switch (iters) { case 1: delta = 75; break; case 2: delta = 60; break; default: delta = 8; if (2*iters + delta < 45) delta = 45 - 2*iters; } below = dlg->filter_aniso_midpoint - delta; above = dlg->filter_aniso_midpoint + delta; contrast_boost_simple(ds_input, below, above, dlg); } } //clear variables delvar: // using = NULL; to get rid of "warning C4701: potentially uninitialized local variable used" compiler message difcoef_lut = NULL; delete [] difcoef_lut; slice_buffer = NULL; delete slice_buffer; //mask_filter->clear(); //mask_filter->data[0][0][0] = 0; //UNREFERENCED_PARAMETER(mask_filter); //to get rid of compiler warning message mask_filter = NULL; delete mask_filter; //create macro files for ImageJ //histogram only <-> (user selects K=1 in GUI) if (dlg->filter_aniso_K != 1) { char file_out_base [512]; sprintf(file_out_base, "%s\\%s(%d_%d_%d).stack.raw", dlg->field_directory_work_str, dlg->field_name_output_str, ds_input->wX, ds_input->wY, ds_input->wZ); //extract macro { char file_out [512]; sprintf(file_out, "%s macro extract.txt", file_out_base); char* outfilemacroextract = file_out; ofstream fout; fout.open(outfilemacroextract, ios::out); // | ios::app); if(!fout.is_open()) { char text[512]; sprintf(text,"histogram: error opening file for saving macro extract: %s", outfilemacroextract); write_to_log(text); return false; } //fill file sprintf(file_out_base, "%s\\%s", dlg->field_directory_work_str, dlg->field_name_output_str); fout << "n = 100;" << endl; fout << "for (i=0; i<=n; i++)" << endl; fout << "{" << endl; fout << " if ( File.exists(\""; fout << file_out_base; fout << "_K" << dlg->filter_aniso_K; fout << "_i#\"+i+\"("; fout << ds_input->wX; fout << "_"; fout << ds_input->wY; fout << "_1).stack.raw\") )"; fout << endl; fout << " {" << endl; fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_K" << dlg->filter_aniso_K; fout << "_i#\"+i+\"("; fout << ds_input->wX; fout << "_"; fout << ds_input->wY; fout << "_1).stack.raw] image=8-bit"; fout << " width="; fout << ds_input->wX; fout << " height="; fout << ds_input->wY; fout << " offset=0 number=1 gap=0\");"; fout << endl; //check if resize needed if ((ds_input->wX > 200+600) && (ds_input->wY > 200+600)) { fout << " makeRectangle(200, 200, 600, 600);" << endl; fout << " run(\"Crop\");" << endl; } // fout << " setColor(255, 255, 255);" << endl; fout << " x=10;" << endl; fout << " y=50;" << endl; fout << " setFont(\"SansSerif\", 40, \"bold\");" << endl; fout << " drawString(\"\"+i+\"\", x, y);" << endl; fout << " }" << endl; fout << "}" << endl; fout << "run(\"Convert Images to Stack\");" << endl; //fout << "saveAs(\"Raw Data\", \""; // fout << file_out_base; // fout << "_crop_600x600.raw\");" << endl; fout << "saveAs(\"AVI... \", \""; fout << file_out_base; fout << "_crop_600x600.avi\");" << endl; fout << endl; fout.close(); // search and replace "\" by "\\" //file buffer char filebuffer [65535]; unsigned int counter; { fstream file; file.open(outfilemacroextract, ios::in | ios::binary ); if(!file.is_open()) { char text[512]; sprintf(text,"histogram: error opening file for saving macro extract: %s", outfilemacroextract); write_to_log(text); return false; } file.seekg (0,ios::beg); counter =0; while (!file.eof()) { filebuffer[counter] = (char)file.get(); counter ++; } file.close(); } { fstream file; file.open(outfilemacroextract, ios::out | ios::binary ); file.seekp (0,ios::beg); unsigned int i; for (i=0; i < counter-1; i++) { file << filebuffer[i]; if (filebuffer[i] == 92) // ASCII(\) = 92 file << "\\"; } file.close(); } } //hist macro sprintf(file_out_base, "%s\\%s(%d_%d_%d).stack.raw", dlg->field_directory_work_str, dlg->field_name_output_str, ds_input->wX, ds_input->wY, ds_input->wZ); { char file_out [512]; sprintf(file_out, "%s macro hist.txt", file_out_base); char* outfilemacrohist = file_out; ofstream fout; fout.open(outfilemacrohist, ios::out); // | ios::app); if(!fout.is_open()) { char text[512]; sprintf(text,"histogram: error opening file for saving macro hist: %s", outfilemacrohist); write_to_log(text); return false; } //fill file fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_hist_260x266.raw] image=8-bit width=260 height=266 offset=0 number=1 gap=0\");"; fout << endl; fout << " setColor(192, 192, 192);" << endl; fout << " x=10; y=40;" << endl; fout << " setFont(\"SansSerif\", 30, \"bold\");" << endl; fout << " drawString(\"0\", x, y);" << endl; fout << endl; //case for variable K (60, 30, then..) if (dlg->filter_aniso_K_var) { //60 fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_K60_i#1_hist_260x266.raw] image=8-bit width=260 height=266 offset=0 number=1 gap=0\");"; fout << endl; fout << " drawString(\"1\", x, y);" << endl; fout << endl; //30 fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_K30_i#2_hist_260x266.raw] image=8-bit width=260 height=266 offset=0 number=1 gap=0\");"; fout << endl; fout << " drawString(\"2\", x, y);" << endl; fout << endl; //the rest fout << "n = 100;" << endl; fout << "for (i=3; i<=n; i++)" << endl; fout << "{" << endl; fout << " if ( File.exists(\""; fout << file_out_base; fout << "_K"; fout << dlg->filter_aniso_K; fout << "_i#\"+i+\"_hist_260x266.raw\") )"; fout << endl; fout << " {" << endl; fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_K"; fout << dlg->filter_aniso_K; fout << "_i#\"+i+\"_hist_260x266.raw] image=8-bit width=260 height=266 offset=0 number=1 gap=0\");"; fout << endl; fout << " drawString(\"\"+i+\"\", x, y);" << endl; fout << " }" << endl; fout << "}" << endl; fout << "run(\"Convert Images to Stack\");" << endl; //fout << "saveAs(\"Raw Data\", \""; // fout << file_out_base; // fout << "_hist_260x266.raw\");" << endl; fout << "saveAs(\"AVI... \", \""; fout << file_out_base; fout << "_hist_260x266.avi\");" << endl; fout << endl; fout.close(); } else //K non var { //the rest fout << "n = 100;" << endl; fout << "for (i=1; i<=n; i++)" << endl; fout << "{" << endl; fout << " if ( File.exists(\""; fout << file_out_base; fout << "_K"; fout << dlg->filter_aniso_K; fout << "_i#\"+i+\"_hist_260x266.raw\") )"; fout << endl; fout << " {" << endl; fout << "run(\"Raw...\", \"open=["; fout << file_out_base; fout << "_K"; fout << dlg->filter_aniso_K; fout << "_i#\"+i+\"_hist_260x266.raw] image=8-bit width=260 height=266 offset=0 number=1 gap=0\");"; fout << endl; fout << " drawString(\"\"+i+\"\", x, y);" << endl; fout << " }" << endl; fout << "}" << endl; fout << "run(\"Convert Images to Stack\");" << endl; //fout << "saveAs(\"Raw Data\", \""; // fout << file_out_base; // fout << "_hist_260x266.raw\");" << endl; fout << "saveAs(\"AVI... \", \""; fout << file_out_base; fout << "_hist_260x266.avi\");" << endl; fout << endl; fout.close(); } // search and replace "\" by "\\" //file buffer char filebuffer [65535]; unsigned int counter; { fstream file; file.open(outfilemacrohist, ios::in | ios::binary ); if(!file.is_open()) { char text[512]; sprintf(text,"histogram: error opening file for saving macro hist: %s", outfilemacrohist); write_to_log(text); return false; } file.seekg (0,ios::beg); counter =0; while (!file.eof()) { filebuffer[counter] = (char)file.get(); counter ++; } file.close(); } { fstream file; file.open(outfilemacrohist, ios::out | ios::binary ); file.seekp (0,ios::beg); unsigned int i; for (i=0; i < counter-1; i++) { file << filebuffer[i]; if (filebuffer[i] == 92) // ASCII(\) = 92 file << "\\"; } file.close(); } } } if (!dlg->stop_request) return true; else return false; } //bool convolve_downsample(DS_cube* ds_input, void * dlg_share) //{ //perform Convolution Filter to downsample cube // PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; // ////determine mask size !!! different than usual mask // unsigned short int mask_size=0; // float radius = dlg->filter_convolve_radius; // mask_size = (short int)(radius) + 1; // // if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) // return false; //sample too small // ////display message in filter window // char win_txt [512]; // sprintf(win_txt, "Convolution Filter - Downsampling by %d...", mask_size); // dlg->txt_field_progress.SetWindowText(win_txt); // // float mask_total = pow(mask_size, 3); // unsigned short int Z2 = (unsigned short int)(ds_input->wZ / mask_size); // //unsigned short int Y2 = (unsigned short int)(ds_input->wY / mask_size); // //unsigned short int X2 = (unsigned short int)(ds_input->wX / mask_size); // ////setup the progress bar parameters // dlg->pbar_simulation_progress.SetRange(0, Z2); // dlg->pbar_simulation_progress.SetPos(0); // dlg->pbar_simulation_progress.SetStep(1); // ////average all slices // { // DS_cube* cube_tot = new DS_cube(1, ds_input->wY, ds_input->wX); // unsigned short int i, j, k; // for (i=0; i < ds_input->wX; i++) // { // sprintf(win_txt, "Averaging (%d / %d)", i, ds_input->wX); // dlg->txt_field_progress.SetWindowText(win_txt); // // for (j=0; j < ds_input->wY; j++) // { // unsigned long long tot = 0; // for (k = 0; k < ds_input->wZ; k++) // if (ds_input->data[k][j][i] < 91) // tot += ds_input->data[k][j][i]; // else // tot += 90; // cube_tot->data[0][j][i] = (unsigned char) ((long double)tot / (long double)ds_input->wZ + .5); // } // } // //save results // char file_out [512]; // sprintf(file_out, "%s_av", dlg->field_directory_and_file_input_str); // cube_tot->write_to_RAW_STACK(file_out); // return false; // } // ////downsample slices // for (unsigned short int depth = 0; depth < Z2; depth++) // { // unsigned short int d2 = mask_size * depth; // //display message in filter window // sprintf(win_txt, "Convolution Filter - Downsampling by %d (depth %d of %d)", mask_size, depth+1, Z2); // dlg->txt_field_progress.SetWindowText(win_txt); // //update progress bar // dlg->pbar_simulation_progress.StepIt(); // // unsigned short int x, y; // for (x=0; x < (unsigned short int)(ds_input->wX / mask_size); x++) // for (y=0; y < (unsigned short int)(ds_input->wY / mask_size); y++) // { // unsigned short int x2 = mask_size * x; // unsigned short int y2 = mask_size * y; // //convolution // float total = 0; // for (unsigned char i=0; i<mask_size; i++) // for (unsigned char j=0; j<mask_size; j++) // for (unsigned char k=0; k<mask_size; k++) // total += ds_input->data[d2+k][y2+j][x2+i]; // //save value // ds_input->data[depth][y][x] = (unsigned char)(total/mask_total+0.5); // } // // if (dlg->stop_request) // { // return false; // } // } ////clear variables //delvar: // return true; //} // bool convolve_edge_detect(DS_cube* ds_input, void * dlg_share) {//perform Convolution Filter for anisotropic_diffusion //display message in filter window PoroMediaDialogShare* dlg = (PoroMediaDialogShare*) dlg_share; char win_txt [512]; sprintf(win_txt, "Convolution Filter - 3D edge detector..."); dlg->txt_field_progress.SetWindowText(win_txt); //determine mask size //float radius = dlg->filter_convolve_radius; float radius = (float)1.9; unsigned short int mask_radius = 1; mask_radius = (short int)(radius); unsigned short int mask_size; mask_size = 2 * mask_radius + 1 ; if ((ds_input->wX < mask_size) || (ds_input->wY < mask_size) || (ds_input->wZ < mask_size)) return false; //sample too small //allocate memory for the mask filter DS_c_cube* sobelx = new DS_c_cube(mask_size, mask_size, mask_size); //DS_si_cube* sobely = new DS_si_cube(mask_size, mask_size, mask_size); //DS_si_cube* sobelz = new DS_si_cube(mask_size, mask_size, mask_size); //laplacian DS_c_cube* lapla = new DS_c_cube(mask_size, mask_size, mask_size); //set values sobelx->data[0][0][0] = 1; sobelx->data[0][0][1] = 0; sobelx->data[0][0][2] = -1; sobelx->data[0][1][0] = 2; sobelx->data[0][1][1] = 0; sobelx->data[0][1][2] = -2; sobelx->data[0][2][0] = -1; sobelx->data[0][2][1] = 0; sobelx->data[0][2][2] = -1; sobelx->data[1][0][0] = 2; sobelx->data[1][0][1] = 0; sobelx->data[1][0][2] = -2; sobelx->data[1][1][0] = 4; sobelx->data[1][1][1] = 0; sobelx->data[1][1][2] = -4; sobelx->data[1][2][0] = 2; sobelx->data[1][2][1] = 0; sobelx->data[1][2][2] = -2; sobelx->data[2][0][0] = 1; sobelx->data[2][0][1] = 0; sobelx->data[2][0][2] = -1; sobelx->data[2][1][0] = 2; sobelx->data[2][1][1] = 0; sobelx->data[2][1][2] = -2; sobelx->data[2][2][0] = -1; sobelx->data[2][2][1] = 0; sobelx->data[2][2][2] = -1; //uniform filling { unsigned short int i,j,k; float radius_sq = radius*radius; for (i=0; i<mask_size; i++) for (j=0; j<mask_size; j++) for (k=0; k<mask_size; k++) if ( pow((float)(i-mask_radius),2) + pow((float)(j-mask_radius),2) + pow((float)(k-mask_radius),2) < radius_sq) lapla->data[k][j][i] = -1; lapla->data[mask_radius][mask_radius][mask_radius] = (char)pow((float)mask_size, 3) - 1 ; } //sobely = rot90x(sobelx); //sobelz = rot90y(sobelx); //allocate mem for cube gradients //store angle values DS_c_cube* gradientx = new DS_c_cube(ds_input->wZ, ds_input->wY, ds_input->wX); DS_c_cube* gradienty = new DS_c_cube(ds_input->wZ, ds_input->wY, ds_input->wX); DS_c_cube* gradientz = new DS_c_cube(ds_input->wZ, ds_input->wY, ds_input->wX); DS_usi_cube* gradientmagnitude = new DS_usi_cube(ds_input->wZ, ds_input->wY, ds_input->wX); DS_cube* gradient_tokeep = new DS_cube(ds_input->wZ, ds_input->wY, ds_input->wX); //could be cDS DS_i_cube* laplacianmagnitude = new DS_i_cube(ds_input->wZ, ds_input->wY, ds_input->wX); //laplacian //set values gradient to zero around the borders of the cube { unsigned short int i,j,k; for (i=0; i<=mask_radius; i++) for (j=0; j<gradient_tokeep->wY; j++) for (k=0; k<gradient_tokeep->wZ; k++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } for (i=gradient_tokeep->wX-mask_radius-1; i < gradient_tokeep->wX; i++) for (j=0; j<gradient_tokeep->wY; j++) for (k=0; k<gradient_tokeep->wZ; k++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } for (j=0; j<=mask_radius-1; j++) for (i=0; i < gradient_tokeep->wX; i++) for (k=0; k<gradient_tokeep->wZ; k++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } for (j=gradient_tokeep->wY-mask_radius-1; j<gradient_tokeep->wY; j++) for (i=0; i < gradient_tokeep->wX; i++) for (k=0; k<gradient_tokeep->wZ; k++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } for (k=0; k<=mask_radius; k++) for (i=0; i < gradient_tokeep->wX; i++) for (j=0; j<gradient_tokeep->wY; j++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } for (k=gradient_tokeep->wZ-mask_radius-1; k<gradient_tokeep->wZ; k++) for (i=0; i < gradient_tokeep->wX; i++) for (j=0; j<gradient_tokeep->wY; j++) { gradient_tokeep->data[k][j][i] = 0; gradientmagnitude->data[k][j][i] = 0; gradientx->data[k][j][i] = 0; gradienty->data[k][j][i] = 0; gradientz->data[k][j][i] = 0; laplacianmagnitude->data[k][j][i] = 0; } } //setup the progress bar parameters dlg->pbar_simulation_progress.SetRange(0, ds_input->wZ - mask_size + 1); dlg->pbar_simulation_progress.SetPos(0); dlg->pbar_simulation_progress.SetStep(1); //calculate gradient results for the first (mask_size - 1) slices for (unsigned short int depth = 0; depth <= ds_input->wZ - mask_size; depth++) { //display message in filter window sprintf(win_txt, "Convolution Filter - 3D edge detector... (depth %d of %d)", depth+1, ds_input->wZ - mask_size + 1); dlg->txt_field_progress.SetWindowText(win_txt); //update progress bar dlg->pbar_simulation_progress.StepIt(); //float sqrt3 = sqrt(3.0); unsigned short int x,y; for (x=0; x <= ds_input->wX - mask_size; x++) for (y=0; y <= ds_input->wY - mask_size; y++) { //convolution float totalx = 0; float totaly = 0; float totalz = 0; float lapl = 0; for (unsigned char i=0; i<mask_size; i++) for (unsigned char j=0; j<mask_size; j++) for (unsigned char k=0; k<mask_size; k++) { totalx += sobelx->data[k][j][i] * ds_input->data[depth+k][y+j][x+i]; totaly += sobelx->data[k][i][j] * ds_input->data[depth+k][y+j][x+i]; totalz += sobelx->data[i][j][k] * ds_input->data[depth+k][y+j][x+i]; lapl += lapla->data[i][j][k] * ds_input->data[depth+k][y+j][x+i]; } //save value float total = sqrt(totalx*totalx + totaly*totaly + totalz*totalz); gradientmagnitude->data[depth+mask_radius][y+mask_radius][x+mask_radius] = (unsigned short int)(total + .5); total *= (float).38; gradientx->data[depth][y+mask_radius][x+mask_radius] = (char)((totalx > total) - (totalx < -total)); gradienty->data[depth][y+mask_radius][x+mask_radius] = (char)((totaly > total) - (totaly < -total)); gradientz->data[depth][y+mask_radius][x+mask_radius] = (char)((totalz > total) - (totalz < -total)); laplacianmagnitude->data[depth][y+mask_radius][x+mask_radius] = (int)lapl; } if (dlg->stop_request) goto delvar; } //non-max suppression (canny type) { //find non-max unsigned short int i,j,k; for (i=mask_radius; i < gradient_tokeep->wX-mask_radius; i++) for (j=mask_radius; j < gradient_tokeep->wY-mask_radius; j++) for (k=mask_radius; k < gradient_tokeep->wZ-mask_radius; k++) { gradient_tokeep->data[k][j][i] = \ ( \ ( gradientmagnitude->data[k][j][i] \ > gradientmagnitude->data[k+gradientz->data[k][j][i]][j+gradienty->data[k][j][i]][i+gradientx->data[k][j][i]] ) \ || ( pow((float)gradientx->data[k][j][i] - gradientx->data[k+gradientz->data[k][j][i]][j+gradienty->data[k][j][i]][i+gradientx->data[k][j][i]], 2) \ + pow((float)gradienty->data[k][j][i] - gradienty->data[k+gradientz->data[k][j][i]][j+gradienty->data[k][j][i]][i+gradientx->data[k][j][i]], 2) \ + pow((float)gradientz->data[k][j][i] - gradientz->data[k+gradientz->data[k][j][i]][j+gradienty->data[k][j][i]][i+gradientx->data[k][j][i]], 2) \ > 2) \ ); gradient_tokeep->data[k][j][i] = ( gradient_tokeep->data[k][j][i] ) && \ ( \ ( gradientmagnitude->data[k][j][i] \ > gradientmagnitude->data[k-gradientz->data[k][j][i]][j-gradienty->data[k][j][i]][i-gradientx->data[k][j][i]] ) \ || ( (pow((float)gradientx->data[k][j][i] - gradientx->data[k-gradientz->data[k][j][i]][j-gradienty->data[k][j][i]][i-gradientx->data[k][j][i]], 2) \ + pow((float)gradienty->data[k][j][i] - gradienty->data[k-gradientz->data[k][j][i]][j-gradienty->data[k][j][i]][i-gradientx->data[k][j][i]], 2) \ + pow((float)gradientz->data[k][j][i]-gradientz->data[k-gradientz->data[k][j][i]][j-gradienty->data[k][j][i]][i-gradientx->data[k][j][i]], 2)) \ > 2) \ ); } //save tokeep char file_out [512]; sprintf(file_out, "%s_tokeep", dlg->field_directory_and_file_input_str); gradient_tokeep->write_to_RAW_STACK(file_out); //update gradientmagnitude for (i=0; i < gradient_tokeep->wX; i++) for (j=0; j < gradient_tokeep->wY; j++) for (k=0; k < gradient_tokeep->wZ; k++) gradientmagnitude->data[k][j][i] = gradientmagnitude->data[k][j][i] * gradient_tokeep->data[k][j][i]; } //save result (gradientmagnitude) { unsigned short int i,j,k; //find max value float max = 0; for (i=mask_radius; i<gradientmagnitude->wX-mask_radius; i++) for (j=mask_radius; j<gradientmagnitude->wY-mask_radius; j++) for (k=mask_radius; k<gradientmagnitude->wZ-mask_radius; k++) if (gradientmagnitude->data[k][j][i] > max) max = (float)gradientmagnitude->data[k][j][i]; max /= (float)255; //rescale & save for (i=0; i<gradientmagnitude->wX; i++) for (j=0; j<gradientmagnitude->wY; j++) for (k=0; k<gradientmagnitude->wZ; k++) ds_input->data[k][j][i] = (unsigned char)( (float)gradientmagnitude->data[k][j][i] / max ); } //clear variables delvar: delete sobelx; delete gradientx; delete gradienty; delete gradientz; delete gradientmagnitude; delete gradient_tokeep; if (dlg->stop_request) return false; else return true; }