From 6f40c95c083464d45205523f6ea75b5197f26a55 Mon Sep 17 00:00:00 2001 From: Gregor Thalhammer Date: Sat, 12 May 2012 21:15:55 +0200 Subject: [PATCH] unwrap: moved global variables (wrap_around) to params struct --- ...wrapper_with_mask_and_wrap_around_option.c | 1418 +++++++++-------- unwrap2D/test_unwrap.py | 9 +- 2 files changed, 724 insertions(+), 703 deletions(-) diff --git a/unwrap2D/Miguel_2D_unwrapper_with_mask_and_wrap_around_option.c b/unwrap2D/Miguel_2D_unwrapper_with_mask_and_wrap_around_option.c index 99dd9cce..d8411c81 100644 --- a/unwrap2D/Miguel_2D_unwrapper_with_mask_and_wrap_around_option.c +++ b/unwrap2D/Miguel_2D_unwrapper_with_mask_and_wrap_around_option.c @@ -1,704 +1,722 @@ -// 2D phase unwrapping, modified for inclusion in scipy by Gregor Thalhammer - -//This program was written by Munther Gdeisat and Miguel Arevallilo HerraŽez to program the two-dimensional unwrapper -//entitled "Fast two-dimensional phase-unwrapping algorithm based on sorting by -//reliability following a noncontinuous path" -//by Miguel Arevallilo HerraŽez, David R. Burton, Michael J. Lalor, and Munther A. Gdeisat -//published in the Journal Applied Optics, Vol. 41, No. 35, pp. 7437, 2002. -//This program was written by Munther Gdeisat, Liverpool John Moores University, United Kingdom. -//Date 26th August 2007 -//The wrapped phase map is assumed to be of floating point data type. The resultant unwrapped phase map is also of floating point type. -//The mask is of byte data type. -//When the mask is 255 this means that the pixel is valid -//When the mask is 0 this means that the pixel is invalid (noisy or corrupted pixel) -//This program takes into consideration the image wrap around problem encountered in MRI imaging. - -//TODO stdlib instead of malloc.h ? (calloc?) -#include -#include -#include - -//TODO: remove global variables -//TODO: make thresholds independent - -static float PI = 3.141592654; -static float TWOPI = 6.283185307; - - -int x_connectivity = 1; -int y_connectivity = 1; -int No_of_edges = 0; - -//PIXELM information -struct pixelm -{ - int increment; //No. of 2*pi to add to the pixel to unwrap it - int number_of_pixels_in_group; //No. of pixel in the pixel group - float value; //value of the pixel - float reliability; - unsigned char input_mask; //0 pixel is masked. 255 pixel is not masked - unsigned char extended_mask; //0 pixel is masked. 255 pixel is not masked - int group; //group No. - int new_group; - struct pixelm *head; //pointer to the first pixel in the group in the linked list - struct pixelm *last; //pointer to the last pixel in the group - struct pixelm *next; //pointer to the next pixel in the group -}; - -typedef struct pixelm PIXELM; - - -//the EDGE is the line that connects two pixels. -//if we have S pixels, then we have S horizontal edges and S vertical edges -struct edge -{ - float reliab; //reliabilty of the edge and it depends on the two pixels +// 2D phase unwrapping, modified for inclusion in scipy by Gregor Thalhammer + +//This program was written by Munther Gdeisat and Miguel Arevallilo HerraŽez to program the two-dimensional unwrapper +//entitled "Fast two-dimensional phase-unwrapping algorithm based on sorting by +//reliability following a noncontinuous path" +//by Miguel Arevallilo HerraŽez, David R. Burton, Michael J. Lalor, and Munther A. Gdeisat +//published in the Journal Applied Optics, Vol. 41, No. 35, pp. 7437, 2002. +//This program was written by Munther Gdeisat, Liverpool John Moores University, United Kingdom. +//Date 26th August 2007 +//The wrapped phase map is assumed to be of floating point data type. The resultant unwrapped phase map is also of floating point type. +//The mask is of byte data type. +//When the mask is 255 this means that the pixel is valid +//When the mask is 0 this means that the pixel is invalid (noisy or corrupted pixel) +//This program takes into consideration the image wrap around problem encountered in MRI imaging. + +//TODO stdlib instead of malloc.h ? (calloc?) +#include +#include +#include + +//TODO: remove global variables +//TODO: make thresholds independent + +static float PI = 3.141592654; +static float TWOPI = 6.283185307; + +typedef struct +{ + float mod; + int x_connectivity; + int y_connectivity; + int no_of_edges; +} params_t; + + +//int x_connectivity = 1; +//int y_connectivity = 1; +int No_of_edges = 0; + +//PIXELM information +struct pixelm +{ + int increment; //No. of 2*pi to add to the pixel to unwrap it + int number_of_pixels_in_group; //No. of pixel in the pixel group + float value; //value of the pixel + float reliability; + unsigned char input_mask; //0 pixel is masked. 255 pixel is not masked + unsigned char extended_mask; //0 pixel is masked. 255 pixel is not masked + int group; //group No. + int new_group; + struct pixelm *head; //pointer to the first pixel in the group in the linked list + struct pixelm *last; //pointer to the last pixel in the group + struct pixelm *next; //pointer to the next pixel in the group +}; + +typedef struct pixelm PIXELM; + +//the EDGE is the line that connects two pixels. +//if we have S pixels, then we have S horizontal edges and S vertical edges +struct edge +{ + float reliab; //reliabilty of the edge and it depends on the two pixels PIXELM *pointer_1; //pointer to the first pixel PIXELM *pointer_2; //pointer to the second pixel - int increment; //No. of 2*pi to add to one of the pixels to unwrap it with respect to the second -}; - -typedef struct edge EDGE; - -//---------------start quicker_sort algorithm -------------------------------- -#define swap(x,y) {EDGE t; t=x; x=y; y=t;} -#define order(x,y) if (x.reliab > y.reliab) swap(x,y) -#define o2(x,y) order(x,y) -#define o3(x,y,z) o2(x,y); o2(x,z); o2(y,z) - -typedef enum {yes, no} yes_no; - -yes_no find_pivot(EDGE *left, EDGE *right, float *pivot_ptr) -{ - EDGE a, b, c, *p; - - a = *left; - b = *(left + (right - left) /2 ); - c = *right; - o3(a,b,c); - - if (a.reliab < b.reliab) - { - *pivot_ptr = b.reliab; - return yes; - } - - if (b.reliab < c.reliab) - { - *pivot_ptr = c.reliab; - return yes; - } - - for (p = left + 1; p <= right; ++p) - { - if (p->reliab != left->reliab) - { - *pivot_ptr = (p->reliab < left->reliab) ? left->reliab : p->reliab; - return yes; - } - return no; - } -} - -EDGE *partition(EDGE *left, EDGE *right, float pivot) -{ - while (left <= right) - { - while (left->reliab < pivot) - ++left; - while (right->reliab >= pivot) - --right; - if (left < right) - { - swap (*left, *right); - ++left; - --right; - } - } - return left; -} - -void quicker_sort(EDGE *left, EDGE *right) -{ - EDGE *p; - float pivot; - - if (find_pivot(left, right, &pivot) == yes) - { - p = partition(left, right, pivot); - quicker_sort(left, p - 1); - quicker_sort(p, right); - } -} - -//--------------end quicker_sort algorithm ----------------------------------- - -//--------------------start initialse pixels ---------------------------------- -//initialse pixels. See the explination of the pixel class above. -//initially every pixel is assumed to belong to a group consisting of only itself -void initialisePIXELs(float *WrappedImage, unsigned char *input_mask, unsigned char *extended_mask, PIXELM *pixel, int image_width, int image_height) -{ - PIXELM *pixel_pointer = pixel; - float *wrapped_image_pointer = WrappedImage; - unsigned char *input_mask_pointer = input_mask; - unsigned char *extended_mask_pointer = extended_mask; - int i, j; - - for (i=0; i < image_height; i++) - { - for (j=0; j < image_width; j++) - { - //pixel_pointer->x = j; - //pixel_pointer->y = i; - pixel_pointer->increment = 0; - pixel_pointer->number_of_pixels_in_group = 1; - pixel_pointer->value = *wrapped_image_pointer; - pixel_pointer->reliability = 9999999 + rand(); - pixel_pointer->input_mask = *input_mask_pointer; - pixel_pointer->extended_mask = *extended_mask_pointer; - pixel_pointer->head = pixel_pointer; - pixel_pointer->last = pixel_pointer; - pixel_pointer->next = NULL; - pixel_pointer->new_group = 0; - pixel_pointer->group = -1; - pixel_pointer++; - wrapped_image_pointer++; - input_mask_pointer++; - extended_mask_pointer++; - } - } -} -//-------------------end initialise pixels ----------- - -//gamma function in the paper -float wrap(float pixel_value) -{ - float wrapped_pixel_value; - if (pixel_value > PI) wrapped_pixel_value = pixel_value - TWOPI; - else if (pixel_value < -PI) wrapped_pixel_value = pixel_value + TWOPI; - else wrapped_pixel_value = pixel_value; - return wrapped_pixel_value; -} - -// pixelL_value is the left pixel, pixelR_value is the right pixel -int find_wrap(float pixelL_value, float pixelR_value) -{ - float difference; - int wrap_value; - difference = pixelL_value - pixelR_value; - - if (difference > PI) wrap_value = -1; - else if (difference < -PI) wrap_value = 1; - else wrap_value = 0; - - return wrap_value; -} - -void extend_mask(unsigned char *input_mask, unsigned char *extended_mask, int image_width, int image_height) -{ - int i,j; - int image_width_plus_one = image_width + 1; - int image_width_minus_one = image_width - 1; - unsigned char *IMP = input_mask + image_width + 1; //input mask pointer - unsigned char *EMP = extended_mask + image_width + 1; //extended mask pointer - - //extend the mask for the image except borders - for (i=1; i < image_height - 1; ++i) - { - for (j=1; j < image_width - 1; ++j) - { - if ( (*IMP) == 255 && (*(IMP + 1) == 255) && (*(IMP - 1) == 255) && - (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && - (*(IMP - image_width_minus_one) == 255) && (*(IMP - image_width_plus_one) == 255) && - (*(IMP + image_width_minus_one) == 255) && (*(IMP + image_width_plus_one) == 255) ) - { - *EMP = 255; - } - ++EMP; - ++IMP; - } - EMP += 2; - IMP += 2; - } - - if (x_connectivity == 1) - { - //extend the mask for the right border of the image - IMP = input_mask + 2 * image_width - 1; - EMP = extended_mask + 2 * image_width -1; - for (i=1; i < image_height - 1; ++ i) - { - if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && - (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && - (*(IMP - image_width - 1) == 255) && (*(IMP - image_width + 1) == 255) && - (*(IMP + image_width - 1) == 255) && (*(IMP - 2 * image_width + 1) == 255) ) - { - *EMP = 255; - } - EMP += image_width; - IMP += image_width; - } - - //extend the mask for the left border of the image - IMP = input_mask + image_width; - EMP = extended_mask + image_width; - for (i=1; i < image_height - 1; ++i) - { - if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && - (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && - (*(IMP - image_width + 1) == 255) && (*(IMP + image_width + 1) == 255) && - (*(IMP + image_width - 1) == 255) && (*(IMP + 2 * image_width - 1) == 255) ) - { - *EMP = 255; - } - EMP += image_width; - IMP += image_width; - } - } - - if (y_connectivity == 1) - { - //extend the mask for the top border of the image - IMP = input_mask + 1; - EMP = extended_mask + 1; - for (i=1; i < image_width - 1; ++i) - { - if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && - (*(IMP + image_width) == 255) && (*(IMP + image_width * (image_height - 1)) == 255) && - (*(IMP + image_width + 1) == 255) && (*(IMP + image_width - 1) == 255) && - (*(IMP + image_width * (image_height - 1) - 1) == 255) && (*(IMP + image_width * (image_height - 1) + 1) == 255) ) - { - *EMP = 255; - } - EMP++; - IMP++; - } - - //extend the mask for the bottom border of the image - IMP = input_mask + image_width * (image_height - 1) + 1; - EMP = extended_mask + image_width * (image_height - 1) + 1; - for (i=1; i < image_width - 1; ++i) - { - if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && - (*(IMP - image_width) == 255) && (*(IMP - image_width - 1) == 255) && (*(IMP - image_width + 1) == 255) && - (*(IMP - image_width * (image_height - 1) ) == 255) && - (*(IMP - image_width * (image_height - 1) - 1) == 255) && - (*(IMP - image_width * (image_height - 1) + 1) == 255) ) - { - *EMP = 255; - } - EMP++; - IMP++; - } - } -} - -void calculate_reliability(float *wrappedImage, PIXELM *pixel, int image_width, int image_height) -{ - int image_width_plus_one = image_width + 1; - int image_width_minus_one = image_width - 1; - PIXELM *pixel_pointer = pixel + image_width_plus_one; - float *WIP = wrappedImage + image_width_plus_one; //WIP is the wrapped image pointer - float H, V, D1, D2; - int i, j; - - for (i = 1; i < image_height -1; ++i) - { - for (j = 1; j < image_width - 1; ++j) - { - if (pixel_pointer->extended_mask == 255) - { - H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); - V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); - D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); - D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); - pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; - } - pixel_pointer++; - WIP++; - } - pixel_pointer += 2; - WIP += 2; - } - - if (x_connectivity == 1) - { - //calculating the reliability for the left border of the image - PIXELM *pixel_pointer = pixel + image_width; - float *WIP = wrappedImage + image_width; - - for (i = 1; i < image_height - 1; ++i) - { - if (pixel_pointer->extended_mask == 255) - { - H = wrap(*(WIP + image_width - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); - V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); - D1 = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); - D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + 2* image_width - 1)); - pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; - } - pixel_pointer += image_width; - WIP += image_width; - } - - //calculating the reliability for the right border of the image - pixel_pointer = pixel + 2 * image_width - 1; - WIP = wrappedImage + 2 * image_width - 1; - - for (i = 1; i < image_height - 1; ++i) - { - if (pixel_pointer->extended_mask == 255) - { - H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP - image_width_minus_one)); - V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); - D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + 1)); - D2 = wrap(*(WIP - 2 * image_width - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); - pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; - } - pixel_pointer += image_width; - WIP += image_width; - } - } - - if (y_connectivity == 1) - { - //calculating the reliability for the top border of the image - PIXELM *pixel_pointer = pixel + 1; - float *WIP = wrappedImage + 1; - - for (i = 1; i < image_width - 1; ++i) - { - if (pixel_pointer->extended_mask == 255) - { - H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); - V = wrap(*(WIP + image_width*(image_height - 1)) - *WIP) - wrap(*WIP - *(WIP + image_width)); - D1 = wrap(*(WIP + image_width*(image_height - 1) - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); - D2 = wrap(*(WIP + image_width*(image_height - 1) + 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); - pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; - } - pixel_pointer++; - WIP++; - } - - //calculating the reliability for the bottom border of the image - pixel_pointer = pixel + (image_height - 1) * image_width + 1; - WIP = wrappedImage + (image_height - 1) * image_width + 1; - - for (i = 1; i < image_width - 1; ++i) - { - if (pixel_pointer->extended_mask == 255) - { - H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); - V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP -(image_height - 1) * (image_width))); - D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) + 1)); - D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) - 1)); - pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; - } - pixel_pointer++; - WIP++; - } - } -} - -//calculate the reliability of the horizontal edges of the image -//it is calculated by adding the reliability of pixel and the relibility of -//its right-hand neighbour -//edge is calculated between a pixel and its next neighbour -void horizontalEDGEs(PIXELM *pixel, EDGE *edge, int image_width, int image_height) -{ - int i, j; - EDGE *edge_pointer = edge; - PIXELM *pixel_pointer = pixel; - - for (i = 0; i < image_height; i++) - { - for (j = 0; j < image_width - 1; j++) - { - if (pixel_pointer->input_mask == 255 && (pixel_pointer + 1)->input_mask == 255) - { - edge_pointer->pointer_1 = pixel_pointer; - edge_pointer->pointer_2 = (pixel_pointer+1); - edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + 1)->reliability; - edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + 1)->value); - edge_pointer++; - No_of_edges++; - } - pixel_pointer++; - } - pixel_pointer++; - } - //construct edges at the right border of the image - if (x_connectivity == 1) - { - pixel_pointer = pixel + image_width - 1; - for (i = 0; i < image_height; i++) - { - if (pixel_pointer->input_mask == 255 && (pixel_pointer - image_width + 1)->input_mask == 255) - { - edge_pointer->pointer_1 = pixel_pointer; - edge_pointer->pointer_2 = (pixel_pointer - image_width + 1); - edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width + 1)->reliability; - edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width + 1)->value); - edge_pointer++; - No_of_edges++; - } - pixel_pointer+=image_width; - } - } -} - -//calculate the reliability of the vertical edges of the image -//it is calculated by adding the reliability of pixel and the relibility of -//its lower neighbour in the image. -void verticalEDGEs(PIXELM *pixel, EDGE *edge, int image_width, int image_height) -{ - int i, j; - PIXELM *pixel_pointer = pixel; - EDGE *edge_pointer = edge + No_of_edges; - - for (i=0; i < image_height - 1; i++) - { - for (j=0; j < image_width; j++) - { - if (pixel_pointer->input_mask == 255 && (pixel_pointer + image_width)->input_mask == 255) - { - edge_pointer->pointer_1 = pixel_pointer; - edge_pointer->pointer_2 = (pixel_pointer + image_width); - edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + image_width)->reliability; - edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + image_width)->value); - edge_pointer++; - No_of_edges++; - } - pixel_pointer++; - } //j loop - } // i loop - - //construct edges that connect at the bottom border of the image - if (y_connectivity == 1) - { - pixel_pointer = pixel + image_width *(image_height - 1); - for (i = 0; i < image_width; i++) - { - if (pixel_pointer->input_mask == 255 && (pixel_pointer - image_width *(image_height - 1))->input_mask == 255) - { - edge_pointer->pointer_1 = pixel_pointer; - edge_pointer->pointer_2 = (pixel_pointer - image_width *(image_height - 1)); - edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width *(image_height - 1))->reliability; - edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width *(image_height - 1))->value); - edge_pointer++; - No_of_edges++; - } - pixel_pointer++; - } - } -} - -//gather the pixels of the image into groups -void gatherPIXELs(EDGE *edge, int image_width, int image_height) -{ - int k; - PIXELM *PIXEL1; - PIXELM *PIXEL2; - PIXELM *group1; - PIXELM *group2; - EDGE *pointer_edge = edge; - int incremento; - - for (k = 0; k < No_of_edges; k++) - { - PIXEL1 = pointer_edge->pointer_1; - PIXEL2 = pointer_edge->pointer_2; - - //PIXELM 1 and PIXELM 2 belong to different groups - //initially each pixel is a group by it self and one pixel can construct a group - //no else or else if to this if - if (PIXEL2->head != PIXEL1->head) - { - //PIXELM 2 is alone in its group - //merge this pixel with PIXELM 1 group and find the number of 2 pi to add - //to or subtract to unwrap it - if ((PIXEL2->next == NULL) && (PIXEL2->head == PIXEL2)) - { - PIXEL1->head->last->next = PIXEL2; - PIXEL1->head->last = PIXEL2; - (PIXEL1->head->number_of_pixels_in_group)++; - PIXEL2->head=PIXEL1->head; - PIXEL2->increment = PIXEL1->increment-pointer_edge->increment; - } - - //PIXELM 1 is alone in its group - //merge this pixel with PIXELM 2 group and find the number of 2 pi to add - //to or subtract to unwrap it - else if ((PIXEL1->next == NULL) && (PIXEL1->head == PIXEL1)) - { - PIXEL2->head->last->next = PIXEL1; - PIXEL2->head->last = PIXEL1; - (PIXEL2->head->number_of_pixels_in_group)++; - PIXEL1->head = PIXEL2->head; - PIXEL1->increment = PIXEL2->increment+pointer_edge->increment; - } - - //PIXELM 1 and PIXELM 2 both have groups - else - { - group1 = PIXEL1->head; - group2 = PIXEL2->head; - //if the no. of pixels in PIXELM 1 group is larger than the no. of pixels - //in PIXELM 2 group. Merge PIXELM 2 group to PIXELM 1 group - //and find the number of wraps between PIXELM 2 group and PIXELM 1 group - //to unwrap PIXELM 2 group with respect to PIXELM 1 group. - //the no. of wraps will be added to PIXELM 2 grop in the future - if (group1->number_of_pixels_in_group > group2->number_of_pixels_in_group) - { - //merge PIXELM 2 with PIXELM 1 group - group1->last->next = group2; - group1->last = group2->last; - group1->number_of_pixels_in_group = group1->number_of_pixels_in_group + group2->number_of_pixels_in_group; - incremento = PIXEL1->increment-pointer_edge->increment - PIXEL2->increment; - //merge the other pixels in PIXELM 2 group to PIXELM 1 group - while (group2 != NULL) - { - group2->head = group1; - group2->increment += incremento; - group2 = group2->next; - } - } - - //if the no. of pixels in PIXELM 2 group is larger than the no. of pixels - //in PIXELM 1 group. Merge PIXELM 1 group to PIXELM 2 group - //and find the number of wraps between PIXELM 2 group and PIXELM 1 group - //to unwrap PIXELM 1 group with respect to PIXELM 2 group. - //the no. of wraps will be added to PIXELM 1 grop in the future - else - { - //merge PIXELM 1 with PIXELM 2 group - group2->last->next = group1; - group2->last = group1->last; - group2->number_of_pixels_in_group = group2->number_of_pixels_in_group + group1->number_of_pixels_in_group; - incremento = PIXEL2->increment + pointer_edge->increment - PIXEL1->increment; - //merge the other pixels in PIXELM 2 group to PIXELM 1 group - while (group1 != NULL) - { - group1->head = group2; - group1->increment += incremento; - group1 = group1->next; - } // while - - } // else - } //else - } //if - pointer_edge++; - } -} - -//unwrap the image -void unwrapImage(PIXELM *pixel, int image_width, int image_height) -{ - int i; - int image_size = image_width * image_height; - PIXELM *pixel_pointer=pixel; - - for (i = 0; i < image_size; i++) - { - pixel_pointer->value += TWOPI * (float)(pixel_pointer->increment); - pixel_pointer++; - } -} - -//set the masked pixels (mask = 0) to the minimum of the unwrapper phase -void maskImage(PIXELM *pixel, unsigned char *input_mask, int image_width, int image_height) -{ - int image_width_plus_one = image_width + 1; - int image_height_plus_one = image_height + 1; - int image_width_minus_one = image_width - 1; - int image_height_minus_one = image_height - 1; - - PIXELM *pointer_pixel = pixel; - unsigned char *IMP = input_mask; //input mask pointer - float min=99999999.; - int i, j; - int image_size = image_width * image_height; - - //find the minimum of the unwrapped phase - for (i = 0; i < image_size; i++) - { - if ((pointer_pixel->value < min) && (*IMP == 255)) - min = pointer_pixel->value; - - pointer_pixel++; - IMP++; - } - - pointer_pixel = pixel; - IMP = input_mask; - - //set the masked pixels to minimum - for (i = 0; i < image_size; i++) - { - if ((*IMP) == 0) - { - pointer_pixel->value = min; - } - pointer_pixel++; - IMP++; - } -} - -//the input to this unwrapper is an array that contains the wrapped phase map. -//copy the image on the buffer passed to this unwrapper to over-write the unwrapped -//phase map on the buffer of the wrapped phase map. -void returnImage(PIXELM *pixel, float *unwrappedImage, int image_width, int image_height) -{ - int i; - int image_size = image_width * image_height; - float *unwrappedImage_pointer = unwrappedImage; - PIXELM *pixel_pointer = pixel; - - for (i=0; i < image_size; i++) - { - *unwrappedImage_pointer = pixel_pointer->value; - pixel_pointer++; - unwrappedImage_pointer++; - } -} - -//the main function of the unwrapper -int unwrap(float* WrappedImage, float* UnwrappedImage, unsigned char* input_mask, - int image_width, int image_height) -{ - unsigned char *extended_mask; - int image_size = image_height * image_width; - int No_of_Edges_initially = 2 * image_width * image_height; - - extended_mask = (unsigned char *) calloc(image_size, sizeof(unsigned char)); - PIXELM *pixel = (PIXELM *) calloc(image_size, sizeof(PIXELM)); - EDGE *edge = (EDGE *) calloc(No_of_Edges_initially, sizeof(EDGE));; - - extend_mask(input_mask, extended_mask, image_width, image_height); - initialisePIXELs(WrappedImage, input_mask, extended_mask, pixel, image_width, image_height); - calculate_reliability(WrappedImage, pixel, image_width, image_height); - horizontalEDGEs(pixel, edge, image_width, image_height); - verticalEDGEs(pixel, edge, image_width, image_height); - - //sort the EDGEs depending on their reiability. The PIXELs with higher - //relibility (small value) first - quicker_sort(edge, edge + No_of_edges - 1); - - //gather PIXELs into groups - gatherPIXELs(edge, image_width, image_height); - - unwrapImage(pixel, image_width, image_height); - maskImage(pixel, input_mask, image_width, image_height); - - //copy the image from PIXELM structure to the unwrapped phase array - //passed to this function - returnImage(pixel, UnwrappedImage, image_width, image_height); - - free(edge); - free(pixel); - free(extended_mask); - - No_of_edges = 0; + int increment; //No. of 2*pi to add to one of the pixels to unwrap it with respect to the second +}; + +typedef struct edge EDGE; + + +//---------------start quicker_sort algorithm -------------------------------- +#define swap(x,y) {EDGE t; t=x; x=y; y=t;} +#define order(x,y) if (x.reliab > y.reliab) swap(x,y) +#define o2(x,y) order(x,y) +#define o3(x,y,z) o2(x,y); o2(x,z); o2(y,z) + +typedef enum {yes, no} yes_no; + +yes_no find_pivot(EDGE *left, EDGE *right, float *pivot_ptr) +{ + EDGE a, b, c, *p; + + a = *left; + b = *(left + (right - left) /2 ); + c = *right; + o3(a,b,c); + + if (a.reliab < b.reliab) + { + *pivot_ptr = b.reliab; + return yes; + } + + if (b.reliab < c.reliab) + { + *pivot_ptr = c.reliab; + return yes; + } + + for (p = left + 1; p <= right; ++p) + { + if (p->reliab != left->reliab) + { + *pivot_ptr = (p->reliab < left->reliab) ? left->reliab : p->reliab; + return yes; + } + return no; + } +} + +EDGE *partition(EDGE *left, EDGE *right, float pivot) +{ + while (left <= right) + { + while (left->reliab < pivot) + ++left; + while (right->reliab >= pivot) + --right; + if (left < right) + { + swap (*left, *right); + ++left; + --right; + } + } + return left; +} + +void quicker_sort(EDGE *left, EDGE *right) +{ + EDGE *p; + float pivot; + + if (find_pivot(left, right, &pivot) == yes) + { + p = partition(left, right, pivot); + quicker_sort(left, p - 1); + quicker_sort(p, right); + } +} + +//--------------end quicker_sort algorithm ----------------------------------- + +//--------------------start initialse pixels ---------------------------------- +//initialse pixels. See the explination of the pixel class above. +//initially every pixel is assumed to belong to a group consisting of only itself +void initialisePIXELs(float *WrappedImage, unsigned char *input_mask, unsigned char *extended_mask, PIXELM *pixel, int image_width, int image_height) +{ + PIXELM *pixel_pointer = pixel; + float *wrapped_image_pointer = WrappedImage; + unsigned char *input_mask_pointer = input_mask; + unsigned char *extended_mask_pointer = extended_mask; + int i, j; + + for (i=0; i < image_height; i++) + { + for (j=0; j < image_width; j++) + { + //pixel_pointer->x = j; + //pixel_pointer->y = i; + pixel_pointer->increment = 0; + pixel_pointer->number_of_pixels_in_group = 1; + pixel_pointer->value = *wrapped_image_pointer; + pixel_pointer->reliability = 9999999 + rand(); + pixel_pointer->input_mask = *input_mask_pointer; + pixel_pointer->extended_mask = *extended_mask_pointer; + pixel_pointer->head = pixel_pointer; + pixel_pointer->last = pixel_pointer; + pixel_pointer->next = NULL; + pixel_pointer->new_group = 0; + pixel_pointer->group = -1; + pixel_pointer++; + wrapped_image_pointer++; + input_mask_pointer++; + extended_mask_pointer++; + } + } +} +//-------------------end initialise pixels ----------- + +//gamma function in the paper +float wrap(float pixel_value) +{ + float wrapped_pixel_value; + if (pixel_value > PI) wrapped_pixel_value = pixel_value - TWOPI; + else if (pixel_value < -PI) wrapped_pixel_value = pixel_value + TWOPI; + else wrapped_pixel_value = pixel_value; + return wrapped_pixel_value; +} + +// pixelL_value is the left pixel, pixelR_value is the right pixel +int find_wrap(float pixelL_value, float pixelR_value) +{ + float difference; + int wrap_value; + difference = pixelL_value - pixelR_value; + + if (difference > PI) wrap_value = -1; + else if (difference < -PI) wrap_value = 1; + else wrap_value = 0; + + return wrap_value; +} + +void extend_mask(unsigned char *input_mask, unsigned char *extended_mask, + int image_width, int image_height, + params_t *params) +{ + int i,j; + int image_width_plus_one = image_width + 1; + int image_width_minus_one = image_width - 1; + unsigned char *IMP = input_mask + image_width + 1; //input mask pointer + unsigned char *EMP = extended_mask + image_width + 1; //extended mask pointer + + //extend the mask for the image except borders + for (i=1; i < image_height - 1; ++i) + { + for (j=1; j < image_width - 1; ++j) + { + if ( (*IMP) == 255 && (*(IMP + 1) == 255) && (*(IMP - 1) == 255) && + (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && + (*(IMP - image_width_minus_one) == 255) && (*(IMP - image_width_plus_one) == 255) && + (*(IMP + image_width_minus_one) == 255) && (*(IMP + image_width_plus_one) == 255) ) + { + *EMP = 255; + } + ++EMP; + ++IMP; + } + EMP += 2; + IMP += 2; + } + + if (params->x_connectivity == 1) + { + //extend the mask for the right border of the image + IMP = input_mask + 2 * image_width - 1; + EMP = extended_mask + 2 * image_width -1; + for (i=1; i < image_height - 1; ++ i) + { + if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && + (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && + (*(IMP - image_width - 1) == 255) && (*(IMP - image_width + 1) == 255) && + (*(IMP + image_width - 1) == 255) && (*(IMP - 2 * image_width + 1) == 255) ) + { + *EMP = 255; + } + EMP += image_width; + IMP += image_width; + } + + //extend the mask for the left border of the image + IMP = input_mask + image_width; + EMP = extended_mask + image_width; + for (i=1; i < image_height - 1; ++i) + { + if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && + (*(IMP + image_width) == 255) && (*(IMP - image_width) == 255) && + (*(IMP - image_width + 1) == 255) && (*(IMP + image_width + 1) == 255) && + (*(IMP + image_width - 1) == 255) && (*(IMP + 2 * image_width - 1) == 255) ) + { + *EMP = 255; + } + EMP += image_width; + IMP += image_width; + } + } + + if (params->y_connectivity == 1) + { + //extend the mask for the top border of the image + IMP = input_mask + 1; + EMP = extended_mask + 1; + for (i=1; i < image_width - 1; ++i) + { + if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && + (*(IMP + image_width) == 255) && (*(IMP + image_width * (image_height - 1)) == 255) && + (*(IMP + image_width + 1) == 255) && (*(IMP + image_width - 1) == 255) && + (*(IMP + image_width * (image_height - 1) - 1) == 255) && (*(IMP + image_width * (image_height - 1) + 1) == 255) ) + { + *EMP = 255; + } + EMP++; + IMP++; + } + + //extend the mask for the bottom border of the image + IMP = input_mask + image_width * (image_height - 1) + 1; + EMP = extended_mask + image_width * (image_height - 1) + 1; + for (i=1; i < image_width - 1; ++i) + { + if ( (*IMP) == 255 && (*(IMP - 1) == 255) && (*(IMP + 1) == 255) && + (*(IMP - image_width) == 255) && (*(IMP - image_width - 1) == 255) && (*(IMP - image_width + 1) == 255) && + (*(IMP - image_width * (image_height - 1) ) == 255) && + (*(IMP - image_width * (image_height - 1) - 1) == 255) && + (*(IMP - image_width * (image_height - 1) + 1) == 255) ) + { + *EMP = 255; + } + EMP++; + IMP++; + } + } +} + +void calculate_reliability(float *wrappedImage, PIXELM *pixel, + int image_width, int image_height, + params_t *params) +{ + int image_width_plus_one = image_width + 1; + int image_width_minus_one = image_width - 1; + PIXELM *pixel_pointer = pixel + image_width_plus_one; + float *WIP = wrappedImage + image_width_plus_one; //WIP is the wrapped image pointer + float H, V, D1, D2; + int i, j; + + for (i = 1; i < image_height -1; ++i) + { + for (j = 1; j < image_width - 1; ++j) + { + if (pixel_pointer->extended_mask == 255) + { + H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); + V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); + D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); + D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); + pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; + } + pixel_pointer++; + WIP++; + } + pixel_pointer += 2; + WIP += 2; + } + + if (params->x_connectivity == 1) + { + //calculating the reliability for the left border of the image + PIXELM *pixel_pointer = pixel + image_width; + float *WIP = wrappedImage + image_width; + + for (i = 1; i < image_height - 1; ++i) + { + if (pixel_pointer->extended_mask == 255) + { + H = wrap(*(WIP + image_width - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); + V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); + D1 = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); + D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + 2* image_width - 1)); + pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; + } + pixel_pointer += image_width; + WIP += image_width; + } + + //calculating the reliability for the right border of the image + pixel_pointer = pixel + 2 * image_width - 1; + WIP = wrappedImage + 2 * image_width - 1; + + for (i = 1; i < image_height - 1; ++i) + { + if (pixel_pointer->extended_mask == 255) + { + H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP - image_width_minus_one)); + V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width)); + D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + 1)); + D2 = wrap(*(WIP - 2 * image_width - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); + pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; + } + pixel_pointer += image_width; + WIP += image_width; + } + } + + if (params->y_connectivity == 1) + { + //calculating the reliability for the top border of the image + PIXELM *pixel_pointer = pixel + 1; + float *WIP = wrappedImage + 1; + + for (i = 1; i < image_width - 1; ++i) + { + if (pixel_pointer->extended_mask == 255) + { + H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); + V = wrap(*(WIP + image_width*(image_height - 1)) - *WIP) - wrap(*WIP - *(WIP + image_width)); + D1 = wrap(*(WIP + image_width*(image_height - 1) - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one)); + D2 = wrap(*(WIP + image_width*(image_height - 1) + 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one)); + pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; + } + pixel_pointer++; + WIP++; + } + + //calculating the reliability for the bottom border of the image + pixel_pointer = pixel + (image_height - 1) * image_width + 1; + WIP = wrappedImage + (image_height - 1) * image_width + 1; + + for (i = 1; i < image_width - 1; ++i) + { + if (pixel_pointer->extended_mask == 255) + { + H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1)); + V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP -(image_height - 1) * (image_width))); + D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) + 1)); + D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) - 1)); + pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2; + } + pixel_pointer++; + WIP++; + } + } +} + +//calculate the reliability of the horizontal edges of the image +//it is calculated by adding the reliability of pixel and the relibility of +//its right-hand neighbour +//edge is calculated between a pixel and its next neighbour +void horizontalEDGEs(PIXELM *pixel, EDGE *edge, + int image_width, int image_height, + params_t *params) +{ + int i, j; + EDGE *edge_pointer = edge; + PIXELM *pixel_pointer = pixel; + + for (i = 0; i < image_height; i++) + { + for (j = 0; j < image_width - 1; j++) + { + if (pixel_pointer->input_mask == 255 && (pixel_pointer + 1)->input_mask == 255) + { + edge_pointer->pointer_1 = pixel_pointer; + edge_pointer->pointer_2 = (pixel_pointer+1); + edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + 1)->reliability; + edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + 1)->value); + edge_pointer++; + No_of_edges++; + } + pixel_pointer++; + } + pixel_pointer++; + } + //construct edges at the right border of the image + if (params->x_connectivity == 1) + { + pixel_pointer = pixel + image_width - 1; + for (i = 0; i < image_height; i++) + { + if (pixel_pointer->input_mask == 255 && (pixel_pointer - image_width + 1)->input_mask == 255) + { + edge_pointer->pointer_1 = pixel_pointer; + edge_pointer->pointer_2 = (pixel_pointer - image_width + 1); + edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width + 1)->reliability; + edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width + 1)->value); + edge_pointer++; + No_of_edges++; + } + pixel_pointer+=image_width; + } + } +} + +//calculate the reliability of the vertical edges of the image +//it is calculated by adding the reliability of pixel and the relibility of +//its lower neighbour in the image. +void verticalEDGEs(PIXELM *pixel, EDGE *edge, + int image_width, int image_height, + params_t *params) +{ + int i, j; + PIXELM *pixel_pointer = pixel; + EDGE *edge_pointer = edge + No_of_edges; + + for (i=0; i < image_height - 1; i++) + { + for (j=0; j < image_width; j++) + { + if (pixel_pointer->input_mask == 255 && (pixel_pointer + image_width)->input_mask == 255) + { + edge_pointer->pointer_1 = pixel_pointer; + edge_pointer->pointer_2 = (pixel_pointer + image_width); + edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + image_width)->reliability; + edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + image_width)->value); + edge_pointer++; + No_of_edges++; + } + pixel_pointer++; + } //j loop + } // i loop + + //construct edges that connect at the bottom border of the image + if (params->y_connectivity == 1) + { + pixel_pointer = pixel + image_width *(image_height - 1); + for (i = 0; i < image_width; i++) + { + if (pixel_pointer->input_mask == 255 && (pixel_pointer - image_width *(image_height - 1))->input_mask == 255) + { + edge_pointer->pointer_1 = pixel_pointer; + edge_pointer->pointer_2 = (pixel_pointer - image_width *(image_height - 1)); + edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width *(image_height - 1))->reliability; + edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width *(image_height - 1))->value); + edge_pointer++; + No_of_edges++; + } + pixel_pointer++; + } + } +} + +//gather the pixels of the image into groups +void gatherPIXELs(EDGE *edge, int image_width, int image_height) +{ + int k; + PIXELM *PIXEL1; + PIXELM *PIXEL2; + PIXELM *group1; + PIXELM *group2; + EDGE *pointer_edge = edge; + int incremento; + + for (k = 0; k < No_of_edges; k++) + { + PIXEL1 = pointer_edge->pointer_1; + PIXEL2 = pointer_edge->pointer_2; + + //PIXELM 1 and PIXELM 2 belong to different groups + //initially each pixel is a group by it self and one pixel can construct a group + //no else or else if to this if + if (PIXEL2->head != PIXEL1->head) + { + //PIXELM 2 is alone in its group + //merge this pixel with PIXELM 1 group and find the number of 2 pi to add + //to or subtract to unwrap it + if ((PIXEL2->next == NULL) && (PIXEL2->head == PIXEL2)) + { + PIXEL1->head->last->next = PIXEL2; + PIXEL1->head->last = PIXEL2; + (PIXEL1->head->number_of_pixels_in_group)++; + PIXEL2->head=PIXEL1->head; + PIXEL2->increment = PIXEL1->increment-pointer_edge->increment; + } + + //PIXELM 1 is alone in its group + //merge this pixel with PIXELM 2 group and find the number of 2 pi to add + //to or subtract to unwrap it + else if ((PIXEL1->next == NULL) && (PIXEL1->head == PIXEL1)) + { + PIXEL2->head->last->next = PIXEL1; + PIXEL2->head->last = PIXEL1; + (PIXEL2->head->number_of_pixels_in_group)++; + PIXEL1->head = PIXEL2->head; + PIXEL1->increment = PIXEL2->increment+pointer_edge->increment; + } + + //PIXELM 1 and PIXELM 2 both have groups + else + { + group1 = PIXEL1->head; + group2 = PIXEL2->head; + //if the no. of pixels in PIXELM 1 group is larger than the no. of pixels + //in PIXELM 2 group. Merge PIXELM 2 group to PIXELM 1 group + //and find the number of wraps between PIXELM 2 group and PIXELM 1 group + //to unwrap PIXELM 2 group with respect to PIXELM 1 group. + //the no. of wraps will be added to PIXELM 2 grop in the future + if (group1->number_of_pixels_in_group > group2->number_of_pixels_in_group) + { + //merge PIXELM 2 with PIXELM 1 group + group1->last->next = group2; + group1->last = group2->last; + group1->number_of_pixels_in_group = group1->number_of_pixels_in_group + group2->number_of_pixels_in_group; + incremento = PIXEL1->increment-pointer_edge->increment - PIXEL2->increment; + //merge the other pixels in PIXELM 2 group to PIXELM 1 group + while (group2 != NULL) + { + group2->head = group1; + group2->increment += incremento; + group2 = group2->next; + } + } + + //if the no. of pixels in PIXELM 2 group is larger than the no. of pixels + //in PIXELM 1 group. Merge PIXELM 1 group to PIXELM 2 group + //and find the number of wraps between PIXELM 2 group and PIXELM 1 group + //to unwrap PIXELM 1 group with respect to PIXELM 2 group. + //the no. of wraps will be added to PIXELM 1 grop in the future + else + { + //merge PIXELM 1 with PIXELM 2 group + group2->last->next = group1; + group2->last = group1->last; + group2->number_of_pixels_in_group = group2->number_of_pixels_in_group + group1->number_of_pixels_in_group; + incremento = PIXEL2->increment + pointer_edge->increment - PIXEL1->increment; + //merge the other pixels in PIXELM 2 group to PIXELM 1 group + while (group1 != NULL) + { + group1->head = group2; + group1->increment += incremento; + group1 = group1->next; + } // while + + } // else + } //else + } //if + pointer_edge++; + } +} + +//unwrap the image +void unwrapImage(PIXELM *pixel, int image_width, int image_height) +{ + int i; + int image_size = image_width * image_height; + PIXELM *pixel_pointer=pixel; + + for (i = 0; i < image_size; i++) + { + pixel_pointer->value += TWOPI * (float)(pixel_pointer->increment); + pixel_pointer++; + } +} + +//set the masked pixels (mask = 0) to the minimum of the unwrapper phase +void maskImage(PIXELM *pixel, unsigned char *input_mask, int image_width, int image_height) +{ + int image_width_plus_one = image_width + 1; + int image_height_plus_one = image_height + 1; + int image_width_minus_one = image_width - 1; + int image_height_minus_one = image_height - 1; + + PIXELM *pointer_pixel = pixel; + unsigned char *IMP = input_mask; //input mask pointer + float min=99999999.; + int i, j; + int image_size = image_width * image_height; + + //find the minimum of the unwrapped phase + for (i = 0; i < image_size; i++) + { + if ((pointer_pixel->value < min) && (*IMP == 255)) + min = pointer_pixel->value; + + pointer_pixel++; + IMP++; + } + + pointer_pixel = pixel; + IMP = input_mask; + + //set the masked pixels to minimum + for (i = 0; i < image_size; i++) + { + if ((*IMP) == 0) + { + pointer_pixel->value = min; + } + pointer_pixel++; + IMP++; + } +} + +//the input to this unwrapper is an array that contains the wrapped +//phase map. copy the image on the buffer passed to this unwrapper to +//over-write the unwrapped phase map on the buffer of the wrapped +//phase map. +void returnImage(PIXELM *pixel, float *unwrappedImage, int image_width, int image_height) +{ + int i; + int image_size = image_width * image_height; + float *unwrappedImage_pointer = unwrappedImage; + PIXELM *pixel_pointer = pixel; + + for (i=0; i < image_size; i++) + { + *unwrappedImage_pointer = pixel_pointer->value; + pixel_pointer++; + unwrappedImage_pointer++; + } +} + +//the main function of the unwrapper +int unwrap(float* WrappedImage, float* UnwrappedImage, unsigned char* input_mask, + int image_width, int image_height) +{ + unsigned char *extended_mask; + int image_size = image_height * image_width; + int No_of_Edges_initially = 2 * image_width * image_height; + params_t params = {TWOPI, 0, 0, 0}; + + extended_mask = (unsigned char *) calloc(image_size, sizeof(unsigned char)); + PIXELM *pixel = (PIXELM *) calloc(image_size, sizeof(PIXELM)); + EDGE *edge = (EDGE *) calloc(No_of_Edges_initially, sizeof(EDGE));; + + extend_mask(input_mask, extended_mask, image_width, image_height, ¶ms); + initialisePIXELs(WrappedImage, input_mask, extended_mask, pixel, image_width, image_height); + calculate_reliability(WrappedImage, pixel, image_width, image_height, ¶ms); + horizontalEDGEs(pixel, edge, image_width, image_height, ¶ms); + verticalEDGEs(pixel, edge, image_width, image_height, ¶ms); + + //sort the EDGEs depending on their reiability. The PIXELs with higher + //relibility (small value) first + quicker_sort(edge, edge + No_of_edges - 1); + + //gather PIXELs into groups + gatherPIXELs(edge, image_width, image_height); + + unwrapImage(pixel, image_width, image_height); + maskImage(pixel, input_mask, image_width, image_height); + + //copy the image from PIXELM structure to the unwrapped phase array + //passed to this function + returnImage(pixel, UnwrappedImage, image_width, image_height); + + free(edge); + free(pixel); + free(extended_mask); + + No_of_edges = 0; return 1; -} +} diff --git a/unwrap2D/test_unwrap.py b/unwrap2D/test_unwrap.py index cad968cf..c7e1aadd 100644 --- a/unwrap2D/test_unwrap.py +++ b/unwrap2D/test_unwrap.py @@ -6,17 +6,17 @@ from numpy import outer, arange, ones, abs, empty, power, indices def test_unwrap2D(): - nx, ny = 64, 64 + nx, ny = 10, 3 x = np.arange(nx) y = np.arange(ny) x.shape = (1,-1) y.shape = (-1,1) - z = np.exp(1j*x*0.5) * np.exp(1j*y*0.1) + z = np.exp(1j*x*0.2*np.pi) * np.exp(1j*y*0.1*np.pi) phi_w = np.angle(z) mask = 255*np.ones((nx, ny), dtype = np.uint8) phi = unwrap2D(phi_w.astype(np.float32), mask) - return phi_w, phi + return phi_w/(np.pi*2), np.asarray(phi)/(np.pi*2) # class test_unwrap(TestCase): @@ -44,4 +44,7 @@ def test_unwrap2D(): if __name__=="__main__": #NumpyTest().run() + import matplotlib.pyplot as plt p,p2 = test_unwrap2D() + plt.imshow(p2) + plt