commit 1a63497f0436458f8a0bcad1b1ff6d7537292c2d Author: Eddy Date: Fri Oct 10 17:36:02 2014 +0400 first commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..44257d3 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +*~ +*.bak +*.bck +*.o +.hg* diff --git a/B-trees.c b/B-trees.c new file mode 100644 index 0000000..aba8516 --- /dev/null +++ b/B-trees.c @@ -0,0 +1,321 @@ +/* + * B-trees.c - my realisation of binary serch trees + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include +#include +#include + +#include "B-trees.h" + +#ifndef Malloc +#define Malloc my_alloc +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} +#endif +#ifndef FREE +#define FREE(arg) do{free(arg); arg = NULL;}while(0) +#endif + +#ifdef EBUG +#ifndef DBG + #define DBG(...) do{fprintf(stderr, __VA_ARGS__);}while(0) +#endif +#else +#ifndef DBG + #define DBG(...) +#endif +#endif + +static inline uint32_t log_2(const uint32_t x) { + uint32_t y; + asm ( "\tbsr %1, %0\n" + : "=r"(y) + : "r" (x) + ); + return y; +} + +// later this function maybe mode difficult +void bt_destroy(BTree *node){ + FREE(node); +} + +// get tree leaf nearest to v +BTree *bt_get(BTree *root, int v){ + if(!root) return NULL; + int D = root->data; + if(D == v) return root; + if(root->data > v) return root->left; + else return root->right; +} + +// find leaf with v or nearest to it (prev) +// prev is nearest parent or NULL for root +BTree *bt_search(BTree *root, int v, BTree **prev){ + if(!root){ + if(prev) *prev = NULL; + return NULL; + } + BTree *p; + do{ + p = root; + root = bt_get(root, v); + }while(root && root->data != v); + if(prev){ + if(p != root) *prev = p; + else *prev = NULL; + } + return root; +} + +BTree *bt_create_leaf(int v){ + BTree *node = Malloc(1, sizeof(BTree)); + if(!node) return NULL; + node->data = v; + return node; +} + +// insert new leaf (or leave it as is, if exists) +// return new node +BTree *bt_insert(BTree **root, int v){ + BTree *closest = NULL, *node = NULL; + if(root && *root){ + node = bt_search(*root, v, &closest); + if(node) return node; // value exists + } + node = bt_create_leaf(v); + if(!node) return NULL; + if(!closest){ // there's no values at all! + if(root) *root = node; // modify root node + return node; + } + // we have a closest node to our data, so create node and insert it: + int D = closest->data; + if(D > v) // there's no left leaf of closest node, add our node there + closest->left = node; + else + closest->right = node; + return node; +} + +// find mostly right element or NULL if no right children +BTree *max_leaf(BTree *root, BTree **prev){ + if(!root) + return NULL; + BTree *ret = root->left; + if(prev) *prev = root; + while(ret->right){ + if(prev) *prev = ret; + ret = ret->right; + } + return ret; +} +// find mostly left element or NULL if no left children +BTree *min_leaf(BTree *root, BTree **prev){ + if(!root) return NULL; + BTree *ret = root->right; + if(prev) *prev = root; + while(ret->left){ + if(prev) *prev = ret; + ret = ret->left; + } + return ret; +} + +void print_tree(BTree *leaf){ + if(!leaf){ + DBG("null"); + return; + } + DBG("(%d: ", leaf->data); + print_tree(leaf->left); + DBG(", "); + print_tree(leaf->right); + DBG(")"); +} + +// remove element, correct root if it is removed +void bt_remove(BTree **root, int v){ + assert(root); + void conn(BTree *p, BTree *n, BTree *node){ + DBG("connect: %d and %d through %d\n", p ? p->data : -1, + n? n->data : -1, node->data); + if(p->left == node) p->left = n; + else p->right = n; + } + void conn_node(BTree *prev, BTree *next, BTree *node){ + if(!prev) *root = next; // this is a root + else conn(prev, next, node); + bt_destroy(node); + } + void node_subst(BTree *nold, BTree *nnew, BTree *pold, BTree *pnew){ + DBG("Substitute: %d by %d\n", nold-> data, nnew->data); + if(!pold) *root = nnew; // this is a root + else conn(pold, nnew, nold); // connect parent of old element to new + conn(pnew, NULL, nnew); // remove node from it's parent + nnew->left = nold->left; + nnew->right = nold->right; + bt_destroy(nold); + } + BTree *node, *prev, *Lmax, *Rmin, *Lmp, *Rmp; + int L, R; + node = bt_search(*root, v, &prev); + if(!node) return; // there's no node v + if(!node->left){ // there's only right child or none at all + conn_node(prev, node->right, node); + return; + } + if(!node->right){ // the same but like mirror + conn_node(prev, node->left, node); + return; + } + // we have situation that element have both children + // 1) find max from left child and min from right + Lmax = max_leaf(node, &Lmp); + L = Lmax->data; + Rmin = min_leaf(node, &Rmp); + R = Rmin->data; + // 2) now L is left max, R - right min + // we select the nearest to v + if(v - L > R - v){ // R is closest value, substitute node by Rmin + node_subst(node, Rmin, prev, Rmp); + }else{ // substitute node by Lmax + node_subst(node, Lmax, prev, Lmp); + } +} + +// return node need rotation +BTree *rotate_ccw(BTree *node, BTree *prev){ + if(!node) return NULL; + #ifdef EBUG + DBG("\n\ntry to rotate near %d", node->data); + if(prev) DBG(", right element: %d", prev->data); + if(!node->right) DBG("\n"); + #endif + if(!node->right) return node->left; // no right child -> return left + BTree *chld = node->right; + DBG(", child: %d\n", chld->data); + if(prev) prev->left = chld; + node->right = chld->left; + chld->left = node; + if(chld->right){DBG("R:%d\n",chld->right->data); return chld;} // need another rotation + else return node; +} + +void tree_to_vine(BTree **root){ + assert(root); assert(*root); + BTree *node = *root, *prev; + while((*root)->right) *root = (*root)->right; // max element will be a new root + while(node && node->right) + node = rotate_ccw(node, NULL); // go to new root + prev = node; node = prev->left; + while(rotate_ccw(node, prev)){ + node = prev->left; + if(!node->right){ + prev = node; + node = node->left; + } + } +} + +// return node->left +BTree *rotate_cw(BTree *node, BTree *parent){ // rotate cw near node->left + if(!node) return NULL; + DBG("\n\ntry to rotate cw near %d\n", node->data); + #ifdef EBUG + if(parent)DBG("\t\tparent: %d\n", parent->data); + #endif + if(!node->left) return NULL; + BTree *chld = node->left, *rght = chld->right; + if(parent) parent->left = chld; + DBG("\t\t, child: %d\n", chld->data); + chld->right = node; + node->left = rght; + return chld; +} + +// make balanced tree +void vine_to_tree(BTree **root){ + uint32_t depth, i; + BTree *node = *root, *prev; + for(depth = 0; node; depth++, node = node->left); + DBG("depth: %u; ", depth); + depth = log_2(depth); + DBG("log: %d\n", depth); + for(i = 0; i < depth; i++){ + DBG("\nIteration %d\n", i); + node = *root; + prev = NULL; + if((*root)->left) *root = (*root)->left; + while(node && node->left){ + node = rotate_cw(node, prev); + if(!node) break; + prev = node; + node = node->left; + } + } +} + +#ifdef STANDALONE +#define INS(X) do{if(!bt_insert(&root, X)) err(1, "alloc error");}while(0) +int main(void){ + BTree *tp, *root = NULL; + int i, ins[] = {7,3,12,1,5,9,20,0,4,6,13,21}; + void search_bt(){ + printf("\nRoot: %d\nTree", root->data); + print_tree(root); printf("\n\n"); + for(i = 0; i < 22; i++){ + tp = bt_search(root, i, NULL); + printf("%d ", i); + if(!tp) printf("not "); + printf("found\n"); + } + } + for(i = 0; i < 12; i++) + INS(ins[i]); + + search_bt(); + printf("now we delete leafs 4, 12 and 7 (old root!)\n"); + bt_remove(&root, 4); bt_remove(&root, 12); bt_remove(&root, 7); + search_bt(); + printf("add items 2, 4, 19\n"); + INS(2); INS(4); INS(19); + search_bt(); + printf("\nnow make a vine\n"); + tree_to_vine(&root); + search_bt(); + printf("\nAnd balance tree\n"); + vine_to_tree(&root); + search_bt(); + return 0; +} +#endif diff --git a/B-trees.h b/B-trees.h new file mode 100644 index 0000000..028c95a --- /dev/null +++ b/B-trees.h @@ -0,0 +1,41 @@ +/* + * B-trees.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __B_TREES_H__ +#define __B_TREES_H__ + +typedef struct btree_node{ + int data; + struct btree_node *left, *right; +} BTree; + +BTree *bt_get(BTree *root, int v); +BTree *bt_search(BTree *root, int v, BTree **prev); +BTree *bt_create_leaf(int v); +BTree *bt_insert(BTree **root, int v); +BTree *max_leaf(BTree *root, BTree **prev); +BTree *min_leaf(BTree *root, BTree **prev); +void bt_remove(BTree **root, int v); + +void print_tree(BTree *leaf); + +#endif // __B_TREES_H__ diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..dd5252e --- /dev/null +++ b/Makefile @@ -0,0 +1,9 @@ +CC = gcc +LDFLAGS = +DEFINES = -DSTANDALONE -DEBUG +CFLAGS = -Wall -Werror + +all : fifo_lifo bidirectional_list B-trees simple_list + +%: %.c + $(CC) $(LDFLAGS) $(CFLAGS) $(DEFINES) $< -o $@ diff --git a/README b/README new file mode 100644 index 0000000..70ee00e --- /dev/null +++ b/README @@ -0,0 +1,13 @@ +This is a set of small usefull utilites & macros + + +bidirectional_list - simple list with operation of searching, inserting, deleting +B-trees - simple but slowly binary search trees with all main operations +fifo_lifo - simple stack-like queues +simple_list - 1-directional list with functions: add element; delete list +usefull_macros - a lot of different macros & functions + * safe memory allocation & freeing + * coloured output on tty & monochromeous on non-tty + * error/warning/debug macros + * MMAP files into memory + diff --git a/SMSD-1.5/Makefile b/SMSD-1.5/Makefile new file mode 100644 index 0000000..19bac1f --- /dev/null +++ b/SMSD-1.5/Makefile @@ -0,0 +1,22 @@ +PROGRAM = client +LDFLAGS = +SRCS = client.c +CC = gcc +DEFINES = -D_XOPEN_SOURCE=501 +CXX = gcc +CFLAGS = -Wall -Werror $(DEFINES) +OBJS = $(SRCS:.c=.o) +all : $(PROGRAM) clean +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +# some addition dependencies +# %.o: %.c +# $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@ +#$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) +# @touch $@ + +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) diff --git a/SMSD-1.5/README b/SMSD-1.5/README new file mode 100644 index 0000000..ebe58ac --- /dev/null +++ b/SMSD-1.5/README @@ -0,0 +1,3 @@ +Simple CLI interface for motorized linear translator Standa 8MT175-150 moving by SMSD-1.5 driver + +Run ./client & press h for help. \ No newline at end of file diff --git a/SMSD-1.5/client.c b/SMSD-1.5/client.c new file mode 100644 index 0000000..11d5c16 --- /dev/null +++ b/SMSD-1.5/client.c @@ -0,0 +1,460 @@ +/* + * client.c - simple terminal client for operationg with + * Standa's 8MT175-150 translator by SMSD-1.5 driver + * + * Hardware operates in microsterpping mode (1/16), + * max current = 1.2A + * voltage = 12V + * "0" of driver connected to end-switch at opposite from motor side + * switch of motor's side connected to "IN1" + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include // tcsetattr +#include // tcsetattr, close, read, write +#include // ioctl +#include // printf, getchar, fopen, perror +#include // exit +#include // read +#include // read +#include // signal +#include // time +#include // memcpy, strcmp etc +#include // int types +#include // gettimeofday + +#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0) + +//double t0; // start time +static int bus_error = 0; // last error of data output +enum{ + NO_ERROR = 0, // normal execution + CODE_ERR, // error of exexuted program code + BUS_ERR, // data transmission error + COMMAND_ERR, // wrong command + CMD_DATA_ERR, // wrong data of command + UNDEFINED_ERR // something else wrong +}; + +FILE *fout = NULL; // file for messages duplicating +char *comdev = "/dev/ttyUSB0"; +int BAUD_RATE = B9600; +struct termio oldtty, tty; // TTY flags +struct termios oldt, newt; // terminal flags +int comfd = -1; // TTY fd + +int erase_ctrlr(); + +/** + * function for different purposes that need to know time intervals + * @return double value: time in seconds + * +double dtime(){ + double t; + struct timeval tv; + gettimeofday(&tv, NULL); + t = tv.tv_sec + ((double)tv.tv_usec)/1e6; + return t; +}*/ + +/** + * Exit & return terminal to old state + * @param ex_stat - status (return code) + */ +void quit(int ex_stat){ + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state + if(comfd > 0){ + erase_ctrlr(); + ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state + close(comfd); + } + if(fout) fclose(fout); + printf("Exit! (%d)\n", ex_stat); + exit(ex_stat); +} + +/** + * Open & setup TTY, terminal + */ +void tty_init(){ + // terminal without echo + tcgetattr(STDIN_FILENO, &oldt); + newt = oldt; + //newt.c_lflag &= ~(ICANON | ECHO); + newt.c_lflag &= ~(ICANON); + if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0) quit(-2); + printf("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + fprintf(stderr,"Can't use port %s\n",comdev); + quit(1); + } + printf(" OK\nGet current settings...\n"); + if(ioctl(comfd,TCGETA,&oldtty) < 0) quit(-1); // Get settings + tty = oldtty; + tty.c_lflag = 0; // ~(ICANON | ECHO | ECHOE | ISIG) + tty.c_oflag = 0; + tty.c_cflag = BAUD_RATE|CS8|CREAD|CLOCAL | PARENB; // 9.6k, 8N1, RW, ignore line ctrl + tty.c_cc[VMIN] = 0; // non-canonical mode + tty.c_cc[VTIME] = 5; + if(ioctl(comfd,TCSETA,&tty) < 0) quit(-1); // set new mode + printf(" OK\n"); +} + +/** + * Read characters from console without echo + * @return char readed + */ +int read_console(char *buf, size_t len){ + int rb = 0; + ssize_t L = 0; + struct timeval tv; + int retval; + fd_set rfds; + tcsetattr(STDIN_FILENO, TCSANOW, &newt); + FD_ZERO(&rfds); + FD_SET(STDIN_FILENO, &rfds); + tv.tv_sec = 0; tv.tv_usec = 1000; + retval = select(1, &rfds, NULL, NULL, &tv); + if(retval){ + if(FD_ISSET(STDIN_FILENO, &rfds)){ + if(len < 2 || !buf) // command works as simple getchar + rb = getchar(); + else{ // read all data from console buffer + if((L = read(STDIN_FILENO, buf, len)) > 0) rb = (int)L; + } + } + } + //tcsetattr(STDIN_FILENO, TCSANOW, &oldt); + return rb; +} + +/** + * getchar() without echo + * wait until at least one character pressed + * @return character readed + */ +int mygetchar(){ // аналог getchar() без необходимости жать Enter + int ret; + do ret = read_console(NULL, 1); + while(ret == 0); + return ret; +} + +/** + * Read data from TTY + * @param buff (o) - buffer for data read + * @param length - buffer len + * @return amount of readed bytes + */ +size_t read_tty(uint8_t *buff, size_t length){ + ssize_t L = 0; + fd_set rfds; + struct timeval tv; + int retval; + FD_ZERO(&rfds); + FD_SET(comfd, &rfds); + tv.tv_sec = 0; tv.tv_usec = 50000; + retval = select(comfd + 1, &rfds, NULL, NULL, &tv); + if(retval < 1) return 0; + if(FD_ISSET(comfd, &rfds)){ + if((L = read(comfd, buff, length)) < 1){ + fprintf(stderr, "ERROR on bus, exit!\n"); + exit(-4); + } + } + return (size_t)L; +} + +/** + * wait for answer from server + * @param sock - socket fd + * @return 0 in case of error or timeout, 1 in case of socket ready + * +int waittoread(int sock){ + fd_set fds; + struct timeval timeout; + int rc; + timeout.tv_sec = 0; + timeout.tv_usec = 1000; + FD_ZERO(&fds); + FD_SET(sock, &fds); + rc = select(sock+1, &fds, NULL, NULL, &timeout); + if(rc < 0){ + perror("select failed"); + return 0; + } + if(rc > 0 && FD_ISSET(sock, &fds)) return 1; + return 0; +} +*/ + +void help(){ + printf("\n\nUse this commands:\n" + "0\tMove to end-switch 0\n" + "1\tMove to end-switch 1\n" + "Lxxx\tMake xxx steps toward zero's end-switch (0 main infinity)\n" + "Rxxx\tMake xxx steps toward end-switch 1 (0 main infinity)\n" + "S\tStop/start motor when program is running\n" + "A\tRun previous command again or stop when running\n" + "E\tErase previous program from controller's memory\n" + "\n" + ); +} + +void write_tty(char *str, int L){ + ssize_t D = write(comfd, str, L); + if(D != L){ + fprintf(stderr, "ERROR on bus, exit!\n"); + exit(-3); + } +} + +#define dup_pr(...) do{printf(__VA_ARGS__); if(fout) fprintf(fout, __VA_ARGS__);}while(0) + +size_t read_ctrl_command(char *buf, size_t L){ // read data from controller to buffer buf + int i, j; + char *ptr = buf; + size_t R; + memset(buf, 0, L); + for(j = 0; j < L; j++, ptr++){ + R = 0; + for(i = 0; i < 10 && !R; i++){ + R = read_tty((uint8_t*)ptr, 1); + } + if(!R){j--; break;} // nothing to read + if(*ptr == '*') // read only one command + break; + if(*ptr < ' '){ // omit spaces & non-characters + j--; ptr--; + } + } + return (size_t) j + 1; +} + +int send_command(char *cmd){ + int L = strlen(cmd); + size_t R = 0; + char ans[256]; + write_tty(cmd, L); + R = read_ctrl_command(ans, 255); + DBG("readed: %s (cmd: %s, R = %zd, L = %d)\n", ans, cmd, R, L); + if(!R || (strncmp(ans, cmd, L) != 0)){ + fprintf(stderr, "Error: controller doesn't respond (answer: %s)\n", ans); + return 0; + } + R = read_ctrl_command(ans, 255); + DBG("readed: %s\n", ans); + if(!R){ // controller is running or error + fprintf(stderr, "Controller doesn't answer!\n"); + return 0; + } + bus_error = NO_ERROR; + //if(strncasecmp(ptr, "HM") + //else + if( strncmp(ans, "E10*", 4) != 0 && + strncmp(ans, "E14*", 4) != 0 && + strncmp(ans, "E12*", 4) != 0){ + fprintf(stderr, "Error in controller: "); + if(strncmp(ans, "E13*", 4) == 0){ + bus_error = CODE_ERR; + fprintf(stderr, "runtime"); + }else if(strncmp(ans, "E15*", 4) == 0){ + bus_error = BUS_ERR; + fprintf(stderr, "databus"); + }else if(strncmp(ans, "E16*", 4) == 0){ + bus_error = COMMAND_ERR; + fprintf(stderr, "command"); + }else if(strncmp(ans, "E19*", 4) == 0){ + bus_error = CMD_DATA_ERR; + fprintf(stderr, "command data"); + }else{ + bus_error = UNDEFINED_ERR; + fprintf(stderr, "undefined (%s)", ans); + } + fprintf(stderr, " error\n"); + return 0; + } + DBG("ALL OK\n"); + return 1; +} +/* +int send5times(char *cmd){ // sends command 'cmd' up to 5 times (if errors), return 0 in case of false + int N, R = 0; + for(N = 0; N < 5 && !R; N++){ + R = send_command(cmd); + } + return R; +}*/ + + +int erase_ctrlr(){ + char *errmsg = "\n\nCan't erase controller's memory: some errors occured!\n\n"; + #define safely_send(x) do{ if(bus_error != NO_ERROR){ \ + fprintf(stderr, errmsg); return 0;} send_command(x); }while(0) + if(!send_command("LD1*")){ // start writing a program + if(bus_error == COMMAND_ERR){ // motor is moving + printf("Found running program, stop it\n"); + if(send_command("ST1*")) + send_command("LD1*"); + }else{ + fprintf(stderr, "Controller doesn't answer: try to press S or E\n"); + return 1; + } + } + safely_send("BG*"); // move address pointer to beginning + safely_send("DS*"); // turn off motor + safely_send("ED*"); // end of program + if(bus_error != NO_ERROR){ + fprintf(stderr, errmsg); + return 0; + } + return 1; +} + +void con_sig(int rb){ + int stepsN = 0, got_command = 0; + char command[256]; + if(rb < 1) return; + if(rb == 'q') quit(0); // q == exit + if(rb == 'L' || rb == 'R'){ + if(!fgets(command, 255, stdin)){ + fprintf(stderr, "You should give amount of steps after commands 'L' and 'R'\n"); + return; + } + stepsN = atoi(command) * 16; // microstepping + if(stepsN < 0 || stepsN > 10000000){ + fprintf(stderr, "\n\nSteps amount should be > -1 and < 625000 (0 means infinity)!\n\n"); + return; + } + } + #define Die_on_error(arg) do{if(!send_command(arg)) goto erase_;}while(0) + if(strchr("LR01", rb)){ // command to execute + got_command = 1; + if(!send_command("LD1*")){ // start writing a program + fprintf(stderr, "Error: previous program is running!\n"); + return; + } + Die_on_error("BG*"); // move address pointer to beginning + Die_on_error("EN*"); // enable power + Die_on_error("SD10000*"); // set speed to max (625 steps per second with 1/16) + } + switch(rb){ + case 'h': + help(); + break; + case '0': + Die_on_error("DL*"); + Die_on_error("HM*"); + break; + case '1': + Die_on_error("DR*"); + Die_on_error("ML*"); + break; + case 'R': + Die_on_error("DR*"); + if(stepsN) + sprintf(command, "MV%d*", stepsN); + else + sprintf(command, "MV*"); + Die_on_error(command); + break; + case 'L': + Die_on_error("DL*"); + if(stepsN) + sprintf(command, "MV%d*", stepsN); + else + sprintf(command, "MV*"); + Die_on_error(command); + break; + case 'S': + Die_on_error("PS1*"); + break; + case 'A': + Die_on_error("ST1*"); + break; + case 'E': + erase_ctrlr(); + break; +/* default: + cmd = (uint8_t) rb; + write(comfd, &cmd, 1);*/ + } + if(got_command){ // there was some command: write ending words + Die_on_error("DS*"); // turn off power from motor at end + Die_on_error("ED*"); // signal about command end + Die_on_error("ST1*");// start program + } + return; +erase_: + if(!erase_ctrlr()) quit(1); +} + +/** + * Get integer value from buffer + * @param buff (i) - buffer with int + * @param len - length of data in buffer (could be 2 or 4) + * @return + * +uint32_t get_int(uint8_t *buff, size_t len){ + int i; + printf("read %zd bytes: ", len); + for(i = 0; i < len; i++) printf("0x%x ", buff[i]); + printf("\n"); + if(len != 2 && len != 4){ + fprintf(stdout, "Bad data length!\n"); + return 0xffffffff; + } + uint32_t data = 0; + uint8_t *i8 = (uint8_t*) &data; + if(len == 2) memcpy(i8, buff, 2); + else memcpy(i8, buff, 4); + return data; +}*/ + +int main(int argc, char *argv[]){ + int rb; + uint8_t buff[128]; + size_t L; + if(argc == 2){ + fout = fopen(argv[1], "a"); + if(!fout){ + perror("Can't open output file"); + exit(-1); + } + setbuf(fout, NULL); + } + tty_init(); + signal(SIGTERM, quit); // kill (-15) + signal(SIGINT, quit); // ctrl+C + signal(SIGQUIT, SIG_IGN); // ctrl+\ . + signal(SIGTSTP, SIG_IGN); // ctrl+Z + setbuf(stdout, NULL); + if(!erase_ctrlr()) quit(1); + //t0 = dtime(); + while(1){ + rb = read_console(NULL, 1); + if(rb > 0) con_sig(rb); + L = read_tty(buff, 127); + if(L){ + buff[L] = 0; + printf("%s", buff); + if(fout) fprintf(fout, "%zd\t%s\n", time(NULL), buff); + } + } +} diff --git a/Trinamic/client/Makefile b/Trinamic/client/Makefile new file mode 100644 index 0000000..19bac1f --- /dev/null +++ b/Trinamic/client/Makefile @@ -0,0 +1,22 @@ +PROGRAM = client +LDFLAGS = +SRCS = client.c +CC = gcc +DEFINES = -D_XOPEN_SOURCE=501 +CXX = gcc +CFLAGS = -Wall -Werror $(DEFINES) +OBJS = $(SRCS:.c=.o) +all : $(PROGRAM) clean +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +# some addition dependencies +# %.o: %.c +# $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@ +#$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) +# @touch $@ + +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) diff --git a/Trinamic/client/client.c b/Trinamic/client/client.c new file mode 100644 index 0000000..8f8f42d --- /dev/null +++ b/Trinamic/client/client.c @@ -0,0 +1,495 @@ +/* + * client.c - simple terminal client + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include // tcsetattr +#include // tcsetattr, close, read, write +#include // ioctl +#include // printf, getchar, fopen, perror +#include // exit +#include // read +#include // read +#include // signal +#include // time +#include // memcpy +#include // int types +#include // gettimeofday + +#define BUFLEN 1024 + +double t0; // start time + +uint32_t motor_speed = 10; // motor speed + +FILE *fout = NULL; // file for messages duplicating +char *comdev = "/dev/ttyUSB0"; +int BAUD_RATE = B115200; +struct termio oldtty, tty; // TTY flags +struct termios oldt, newt; // terminal flags +int comfd; // TTY fd + +#define DBG(...) do{fprintf(stderr, __VA_ARGS__);}while(0) + +/** + * function for different purposes that need to know time intervals + * @return double value: time in seconds + * +double dtime(){ + double t; + struct timeval tv; + gettimeofday(&tv, NULL); + t = tv.tv_sec + ((double)tv.tv_usec)/1e6; + return t; +}*/ + +/** + * Exit & return terminal to old state + * @param ex_stat - status (return code) + */ +void quit(int ex_stat){ + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state + ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state + close(comfd); + if(fout) fclose(fout); + printf("Exit! (%d)\n", ex_stat); + exit(ex_stat); +} + +int send_command(uint8_t *ninebytes){ + uint8_t crc = 0; + int i; + printf("send: "); + for(i = 0; i < 8; crc += ninebytes[i++]) + printf("%u,", ninebytes[i]); + ninebytes[8] = crc; + printf("%u\n",ninebytes[8]); + if(9 != write(comfd, ninebytes, 9)){ + perror("Can't write to Trinamic"); + return 0; + } + return 1; +} + + +/** + * Open & setup TTY, terminal + */ +void tty_init(){ + // terminal without echo + tcgetattr(STDIN_FILENO, &oldt); + newt = oldt; + newt.c_lflag &= ~(ICANON | ECHO); + if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0) quit(-2); + printf("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + fprintf(stderr,"Can't use port %s\n",comdev); + quit(1); + } + printf(" OK\nGet current settings...\n"); + if(ioctl(comfd,TCGETA,&oldtty) < 0) quit(-1); // Get settings + tty = oldtty; + tty.c_lflag = 0; // ~(ICANON | ECHO | ECHOE | ISIG) + tty.c_oflag = 0; + tty.c_cflag = BAUD_RATE|CS8|CREAD|CLOCAL; // 9.6k, 8N1, RW, ignore line ctrl + tty.c_cc[VMIN] = 0; // non-canonical mode + tty.c_cc[VTIME] = 5; + if(ioctl(comfd,TCSETA,&tty) < 0) quit(-1); // set new mode + printf(" OK\n"); + tcsetattr(STDIN_FILENO, TCSANOW, &newt); +} + +/** + * getchar() without echo + * wait until at least one character pressed + * @return character readed + * +int mygetchar(){ + int ret; + do ret = read_console(); + while(ret == 0); + return ret; +}*/ + +/** + * read both tty & console + * @param buff (o) - buffer for messages readed from tty or NULL (to omit readed info) + * @param length (io) - buff's length (return readed len or 0) + * @param rb (io) - byte readed from console, -1 if nothing read, NULL to omit console reading + * @return 1 if something was readed here or there + */ +int read_tty_and_console(uint8_t *buff, size_t *length, int *rb){ + ssize_t L; + ssize_t l; + uint8_t buff1[BUFLEN+1]; + size_t buffsz; + struct timeval tv; + int sel, retval = 0; + fd_set rfds; + FD_ZERO(&rfds); + FD_SET(STDIN_FILENO, &rfds); + FD_SET(comfd, &rfds); + tv.tv_sec = 0; tv.tv_usec = 10000; + sel = select(comfd + 1, &rfds, NULL, NULL, &tv); + if(!buff){ + buffsz = BUFLEN; + buff = buff1; + }else{ + buffsz = *length; + } + if(sel > 0){ + if(rb){ + if(FD_ISSET(STDIN_FILENO, &rfds)){ + *rb = getchar(); + retval = 1; + }else{ + *rb = -1; + } + } + if(FD_ISSET(comfd, &rfds)){ + if((L = read(comfd, buff, buffsz)) < 1){ // disconnect or other troubles + fprintf(stderr, "USB error or disconnected!\n"); + quit(1); + }else{ + if(L == 0){ // USB disconnected + fprintf(stderr, "USB disconnected!\n"); + quit(1); + } + + // all OK continue reading + //DBG("readed %zd bytes, try more.. ", L); + buffsz -= L; + usleep(10000); + while(buffsz > 0 && (l = read(comfd, buff+L, buffsz)) > 0){ + L += l; + buffsz -= l; + usleep(10000); + } + //DBG("full len: %zd\n", L); + + if(length) *length = (size_t) L; + retval = 1; + } + }else{ + if(length) *length = 0; + } + } + return retval; +} + +int32_t get_integer(uint8_t *buff){ + int32_t val; + int ii; + uint8_t *valptr = (uint8_t*) &val + 3; + for(ii = 4; ii < 8; ii++) + *valptr-- = buff[ii]; + return val; +} + +void help(){ + printf("Use this commands:\n" + "h\tShow this help\n" + "+/-\tIncrease/decrease speed by 10%% or 1\n" + "M/m\tMove to relative position\n" + "R/L\tRotate motor 1 right/left\n" + "r/l\tRotate motor 2 right/left\n" + "S/s\tStop 1st or 2nd motor\n" + "W/w\tSet ustepping to 1/16\n" + "E/e\tSet ustepping to 1/4\n" + "Z/z\tSet to zero current position\n" + "G/g\tGet current position\n" + "A/a\tGet ALL information about\n" + "q\tQuit\n" + ); +} + +typedef struct{ + char *name; // command name + uint8_t cmd; // command number +} custom_command; + +void get_information(uint8_t *com){ + int i; + uint8_t buff[9]; + size_t L; + custom_command info[] = { + {"target position", 0}, + {"actual position", 1}, + {"target speed", 2}, + {"actual speed", 3}, + {"max positioning speed", 4}, + {"max acceleration", 5}, + {"abs. max current", 6}, + {"standby current", 7}, + {"target pos. reached", 8}, + {"ref. switch status", 9}, + {"right limit switch status", 10}, + {"left limit switch status", 11}, + {"right limit switch disable", 12}, + {"left limit switch disable", 13}, + {"minimum speed", 130}, + {"actual acceleration", 135}, + {"ramp mode", 138}, + {"microstep resolution", 140}, + {"soft stop flag", 149}, + {"ramp divisor", 153}, + {"pulse divisor", 154}, + {"referencing mode", 193}, + {"referencing search speed", 194}, + {"referencing switch speed", 195}, + {"distance end switches", 196}, + {"mixed decay threshold", 203}, + {"freewheeling", 204}, + {"stall detection threshold", 205}, + {"actual load value", 206}, + {"driver error flags", 208}, + {"fullstep threshold", 211}, + {"power down delay", 214}, + {NULL, 0} + }; + printf("\n"); + for(i = 0; info[i].name; i++){ + printf("%s (%u) ", info[i].name, info[i].cmd); + com[2] = info[i].cmd; + if(send_command(com)){ // command send OK + L = 9; + if(read_tty_and_console(buff, &L, NULL)){ // get answer + if(L == 9 && buff[2] > 6){ // answer's length is right && no error in status + int32_t val = get_integer(buff); + printf("= %d\n", val); + //for(ii = 0; ii < 9; ii++) printf("%u,",buff[ii]); + }else printf("L=%zd (stat: %u)\n", L, buff[2]); + }} + printf("\n"); + } +} + +void con_sig(int rb){ + if(rb < 1) return; + uint32_t value = motor_speed; // here is the value to send + /* + * trinamic binary format (bytes): + * 0 - address + * 1 - command + * 2 - type + * 3 - motor/bank number + * 4 \ + * 5 | value (msb first!) + * 6 | + * 7 / + * 8 - checksum (sum of bytes 0..8) + */ + uint8_t command[9] = {1,0,0,0,0,0,0,0,0}; // command to sent + int i, V; + uint8_t *cmd = &command[1]; // byte 1 of message is command itself + if(rb < 1) return; + if(rb == 'q' || rb == 'Q') quit(0); + if(rb == 'h' || rb == 'H'){ + help(); + return; + } + void change_motor_speed(int dir){ + uint32_t Delta = motor_speed / 10, newspd; + int i; + if(Delta == 0) Delta = 1; + if(dir > 0){ + newspd = motor_speed + Delta; + if(newspd > motor_speed) motor_speed = newspd; + else motor_speed = (uint32_t)0xffffffffUL; + }else{ + newspd = motor_speed - Delta; + if(newspd < motor_speed && newspd) motor_speed = newspd; + else motor_speed = 1; + } + printf("\nspeed changed to %u\n", motor_speed); + *cmd = 5; // SAP + command[2] = 4; // max pos speed + uint8_t *spd = (uint8_t*) &motor_speed + 3; + for(i = 4; i < 8; i++) + command[i] = *spd--; + send_command(command); // store new speed + read_tty_and_console(NULL, NULL, NULL); + command[3] = 2; // and for second motor + send_command(command); + read_tty_and_console(NULL, NULL, NULL); + } + if(rb == '+'){ + change_motor_speed(1); + }else if(rb == '-'){ + change_motor_speed(-1); + } + /* + * Now we must analize commands: all letters exept Q/q & H/h are for motors commands + */ + if(rb > 'a'-1 && rb < 'z'+1){ // small letter -> change motor number to 2 + command[3] = 2; + }else if(rb > 'A'-1 && rb < 'Z'+1){ // capital letter -> change to small + rb += 'a' - 'A'; + }else return; // wrong command + switch(rb){ + case 'r': // move motor right + *cmd = 1; // ROR + break; + case 'l': // move motor left + *cmd = 2; // ROL + break; + case 's': // stop motor + *cmd = 3; // STP + break; + case 'w': // 1/16 + *cmd = 5; // SAP + command[2] = 140; // ustep resolution + value = 4; + break; + case 'e': // 1/4 + *cmd = 5; // SAP + command[2] = 140; // ustep resolution + value = 2; + break; + case 'g': // get current position + *cmd = 6; // GAP + command[2] = 1; // actual position + break; + case 'a': // get everything! + *cmd = 6; // GAP + get_information(command); + return; + break; + case 'z': // set to zero current position + *cmd = 5; // SAP + command[2] = 1; // actual position + value = 0; + break; + case 'm': // move to relative position + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to original state + if(scanf("%d", &V) > 0 && V){ + *cmd = 4; // MVP + command[2] = 1; // REL + value = (uint32_t) V; + //if(V > 0) value = (uint32_t) V; + //else value = + } + tcsetattr(STDIN_FILENO, TCSANOW, &newt); // omit echo + break; + } + if(*cmd){ + // copy param to buffer + uint8_t *spd = (uint8_t*) &value + 3; + for(i = 4; i < 8; i++) + command[i] = *spd--; + send_command(command); + } +} + +/** + * Get integer value from buffer + * @param buff (i) - buffer with int + * @param len - length of data in buffer (could be 2 or 4) + * @return + * +uint32_t get_int(uint8_t *buff, size_t len){ + if(len != 2 && len != 4){ + fprintf(stdout, "Bad data length!\n"); + return 0xffffffff; + } + uint32_t data = 0; + uint8_t *i8 = (uint8_t*) &data; + if(len == 2) memcpy(i8, buff, 2); + else memcpy(i8, buff, 4); + return data; +}*/ + +/** + * Copy line by line buffer buff to file removing cmd starting from newline + * @param buffer - data to put into file + * @param cmd - symbol to remove from line startint (if found, change *cmd to (-1) + * or NULL, (-1) if no command to remove + */ +void copy_buf_to_file(uint8_t *buffer, int *cmd){ + uint8_t *buff, *line, *ptr; + if(!cmd || *cmd < 0){ + fprintf(fout, "%s", buffer); + return; + } + buff = (uint8_t*)strdup((char*)buffer), ptr = buff; + do{ + if(!*ptr) break; + if(ptr[0] == (char)*cmd){ + *cmd = -1; + ptr++; + if(ptr[0] == '\n') ptr++; + if(!*ptr) break; + } + line = ptr; + ptr = (uint8_t*)strchr((char*)buff, '\n'); + if(ptr){ + *ptr++ = 0; + //fprintf(fout, "%s\n", line); + }//else + //fprintf(fout, "%s", line); // no newline found in buffer + fprintf(fout, "%s\n", line); + }while(ptr); + free(buff); +} + + +int main(int argc, char *argv[]){ + int rb, oldcmd = -1; + uint8_t buff[BUFLEN+1]; + size_t L; + if(argc == 2){ + fout = fopen(argv[1], "a"); + if(!fout){ + perror("Can't open output file"); + exit(-1); + } + setbuf(fout, NULL); + } + tty_init(); + signal(SIGTERM, quit); // kill (-15) + signal(SIGINT, quit); // ctrl+C + signal(SIGQUIT, SIG_IGN); // ctrl+\ . + signal(SIGTSTP, SIG_IGN); // ctrl+Z + setbuf(stdout, NULL); + //t0 = dtime(); + while(1){ + L = BUFLEN; + if(read_tty_and_console(buff, &L, &rb)){ + if(rb > 0){ + con_sig(rb); + oldcmd = rb; + } + if(L > 0){ + uint8_t *ptr = buff; + int ii; + int32_t val = get_integer(buff); + printf("value = %d (full answer: ", val); + for(ii = L - 1; ii > 0; ii--){ + uint8_t C = *ptr++; + printf("%u", C); + //if(C > 31) printf("(%c)", C); + printf(","); + } + printf("%u)\n", *ptr); + if(fout){ + copy_buf_to_file(buff, &oldcmd); + } + } + } + } +} diff --git a/Trinamic/websock_simple/Makefile b/Trinamic/websock_simple/Makefile new file mode 100644 index 0000000..2704c87 --- /dev/null +++ b/Trinamic/websock_simple/Makefile @@ -0,0 +1,22 @@ +PROGRAM = websocktest +LDFLAGS = $(shell pkg-config --libs libwebsockets) +SRCS = test.c +CC = gcc +DEFINES = -D_XOPEN_SOURCE=501 -DCUR_PATH=\"$(shell pwd)\" +CXX = gcc +CFLAGS = -Wall -Werror -Wextra $(DEFINES) $(shell pkg-config --cflags libwebsockets) +OBJS = $(SRCS:.c=.o) +all : $(PROGRAM) clean +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +# some addition dependencies +# %.o: %.c +# $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@ +#$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) +# @touch $@ + +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) diff --git a/Trinamic/websock_simple/test.c b/Trinamic/websock_simple/test.c new file mode 100644 index 0000000..4378dcb --- /dev/null +++ b/Trinamic/websock_simple/test.c @@ -0,0 +1,183 @@ +#include +#include +#include +#include +#include + +#include + +#include + +#include // tcsetattr +#include // int types +#include // read +#include // read + + +#define _U_ __attribute__((__unused__)) + +int force_exit = 0; + +uint8_t buf[9]; + +char *comdev = "/dev/ttyUSB0"; +int BAUD_RATE = B115200; +int comfd = -1; // TTY fd +int send_command(uint8_t *ninebytes){ + uint8_t crc = 0; + int i; + printf("send: "); + for(i = 0; i < 8; crc += ninebytes[i++]) + printf("%u,", ninebytes[i]); + ninebytes[8] = crc; + printf("%u\n",ninebytes[8]); + if(9 != write(comfd, ninebytes, 9)){ + perror("Can't write to Trinamic"); + return 0; + } + return 1; +} + +void fillbuf(char *command){ + memset(buf, 0, 9); + buf[0] = 1; // controller # + if(command[1] == 'X') buf[3] = 2; // X motor -> #2 + else buf[3] = 0; // Y motor -> #0 + if(command[0] == 'D'){ // start moving + if(command[2] == '+') buf[1] = 1; // ROR + else buf[1] = 2; // ROL + }else{ // stop + buf[1] = 3; // STP + } + buf[7] = 50; // speed = 50 + send_command(buf); +} + +void move_motor(char *where){ + printf("moving %s\n", where); + +} +void stop_motor(char *where){ + printf("stoping %s\n", where); +} + +static int +my_protocol_callback(_U_ struct libwebsocket_context *context, + _U_ struct libwebsocket *wsi, + enum libwebsocket_callback_reasons reason, + _U_ void *user, void *in, _U_ size_t len) +{ + //int n, m; + //unsigned char buf[LWS_SEND_BUFFER_PRE_PADDING + 512 + + // LWS_SEND_BUFFER_POST_PADDING]; + //unsigned char *p = &buf[LWS_SEND_BUFFER_PRE_PADDING]; + + char *msg = (char*) in; + + switch (reason) { + + case LWS_CALLBACK_ESTABLISHED: + printf("New Connection\n"); + break; + + case LWS_CALLBACK_SERVER_WRITEABLE: + break; + + case LWS_CALLBACK_RECEIVE: + fillbuf(msg); + break; + + case LWS_CALLBACK_FILTER_PROTOCOL_CONNECTION: + break; + + default: + break; + } + + return 0; +} +//**************************************************************************// +/* list of supported protocols and callbacks */ +//**************************************************************************// +static struct libwebsocket_protocols protocols[] = { + /* first protocol must always be HTTP handler */ + + { + "XY-protocol", // name + my_protocol_callback, // callback + 30, // per_session_data_size + 10, // max frame size / rx buffer + 0, NULL, 0 + }, + { NULL, NULL, 0, 0, 0, NULL, 0} /* terminator */ +}; +//**************************************************************************// +void sighandler(_U_ int sig){ + close(comfd); + force_exit = 1; +} +//**************************************************************************// +//**************************************************************************// +int main(_U_ int argc, _U_ char **argv) +{ + int n = 0; + struct libwebsocket_context *context; + int opts = 0; + const char *iface = NULL; + int syslog_options = LOG_PID | LOG_PERROR; + //unsigned int oldus = 0; + struct lws_context_creation_info info; + + int debug_level = 7; + + memset(&info, 0, sizeof info); + info.port = 9999; + + signal(SIGINT, sighandler); + + printf("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + fprintf(stderr,"Can't use port %s\n",comdev); + return 1; + } + + + /* we will only try to log things according to our debug_level */ + setlogmask(LOG_UPTO (LOG_DEBUG)); + openlog("lwsts", syslog_options, LOG_DAEMON); + + /* tell the library what debug level to emit and to send it to syslog */ + lws_set_log_level(debug_level, lwsl_emit_syslog); + + info.iface = iface; + info.protocols = protocols; + info.extensions = libwebsocket_get_internal_extensions(); + info.ssl_cert_filepath = NULL; + info.ssl_private_key_filepath = NULL; + + info.gid = -1; + info.uid = -1; + info.options = opts; + + context = libwebsocket_create_context(&info); + if (context == NULL) { + lwsl_err("libwebsocket init failed\n"); + return -1; + } + + n = 0; + while (n >= 0 && !force_exit) { + + n = libwebsocket_service(context, 50); + + };//while n>=0 + + + libwebsocket_context_destroy(context); + + lwsl_notice("libwebsockets-test-server exited cleanly\n"); + + closelog(); + + return 0; +}//main diff --git a/Trinamic/websock_simple/test.html b/Trinamic/websock_simple/test.html new file mode 100644 index 0000000..9bd1449 --- /dev/null +++ b/Trinamic/websock_simple/test.html @@ -0,0 +1,106 @@ + + + + A test + + + +
+ + + + + + + + + + + + +
+

+
No connection
+
+
+ + + diff --git a/Zernike/Makefile b/Zernike/Makefile new file mode 100644 index 0000000..ab301ce --- /dev/null +++ b/Zernike/Makefile @@ -0,0 +1,23 @@ +LOADLIBES = -lm -lgsl -lgslcblas +SRCS = zernike.c zernikeR.c zernike_annular.c +CC = gcc +DEFINES = -D_GNU_SOURCE +#-D_XOPEN_SOURCE=501 +CXX = gcc +CFLAGS = -Wall -Werror $(DEFINES) +OBJS = $(SRCS:.c=.o) +all : zernike btatest +zernike : $(OBJS) test.o + $(CC) $(CFLAGS) test.o $(OBJS) $(LOADLIBES) -o zernike +btatest : $(OBJS) Z-BTA_test.o simple_list.o + $(CC) $(CFLAGS) Z-BTA_test.o simple_list.o $(OBJS) $(LOADLIBES) -o btatest +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(SRCS) + +### +# name1.o : header1.h header2.h ... +test.o zernike.o zernikeR.o zernike_annular.o Z-BTA_test.o : zernike.h +zernike.o zernikeR.o zernike_annular.o : zern_private.h +simple_list.o Z-BTA_test.o : simple_list.h diff --git a/Zernike/README b/Zernike/README new file mode 100644 index 0000000..7655e22 --- /dev/null +++ b/Zernike/README @@ -0,0 +1,26 @@ +Here is my realisation of wavefront decomposition and restoration based on +Zernike polynomials (up to 100th order). + +The base functions are: + +double *zernfunR(int n, int m, int W, int H, double *norm); + calculates Zernike polynomial of order n and angular order m on + equidistance rectangular grid WxH (norm is NULL or value to store + sum(Z_i^2) for normalisation on decomposition. + +double *zernfunNR(int p, int W, int H, double *norm); + the same but in Noll notation + +double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx); + decompose image by Zernike polynomials with order <= Nmax, + Zsz is size of returned array (in Noll notation) + lastidx is last non-zero coefficient + +double *Zcompose(int Zsz, double *Zidxs, int W, int H); + build image based on Zernike coeffs Zidxs + Zsz - size of array in Noll notation + W, H - size of future image + +void set_prec(double val); + set precision of Zernike transforms + diff --git a/Zernike/Z-BTA_test.c b/Zernike/Z-BTA_test.c new file mode 100644 index 0000000..a3bc18d --- /dev/null +++ b/Zernike/Z-BTA_test.c @@ -0,0 +1,402 @@ +/* + * Z-BTA_test.c - simple test for models of hartmannograms + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include +#include +#include +#include +#include +#include // open/close/etc +#include +#include +#include // mmap + +#include "zernike.h" +#include "simple_list.h" + +#define _(...) __VA_ARGS__ +extern char *__progname; + +#define ERR(...) err(1, __VA_ARGS__) +#define ERRX(...) errx(1, __VA_ARGS__) + +// x0 = (*spot)->c.x - AxisX; +// y0 = AxisY - (*spot)->c.y; +double AxisX = 1523., AxisY = 1533.; // BUG from fitsview : xreal = x + AxisX, yreal = y - AxisY + +char *name0 = NULL, *name1 = NULL; // prefocal/postfocal filenames +double distance = -1.; // distance between images +double pixsize = 30.; // pixel size +double ImHeight = 500.; // half of image height in pixels + +// spots for spotlist +typedef struct{ + int id; // spot identificator + double x; // coordinates of center + double y; +} spot; + +// hartmannogram +typedef struct{ + char *filename; // spots-filename + int len; // amount of spots + List *spots; // spotlist +} hartmann; + + +/** + * print usage with optional message & exit with error(1) + * @param fmt ... - printf-like message, no tailing "\n" required! + */ +void usage(char *fmt, ...){ + va_list ap; + va_start(ap, fmt); + printf("\n"); + if (fmt != NULL){ + vprintf(fmt, ap); + printf("\n\n"); + } + va_end(ap); + // "Использование:\t%s опции\n" + printf(_("Usage:\t%s options\n"), + __progname); + // "Опции:\n" + printf(_("Required options:\n")); + printf("\t--prefocal, -0 \t\tfilename with spotslist in " RED "prefocal" OLDCOLOR " image\n"); + printf("\t--postfocal,-1 \t\tfilename with spotslist in " RED "postfocal" OLDCOLOR " image\n"); + printf("\t--distance, -D \tdistance between images in millimeters\n"); + printf(_("Unnesessary options:\n")); + printf("\t--pixsize, -p \tpixel size in microns (default: 30)\n"); + exit(1); +} + +/** + * converts string into double number + * @param num (o) - result of conversion + * @param str (i) - string with number + * @return num of NULL in case of error + */ +double *myatod(double *num, const char *str){ + double tmp; + char *endptr; + if(!num) return NULL; + if(!str) return NULL; + tmp = strtod(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0') + return NULL; + *num = tmp; + return num; +} + +/** + * Parce arguments of main() & fill global options + */ +void parse_args(int argc, char **argv){ + int i; + char short_options[] = "0:1:D:p:"; // all short equivalents + struct option long_options[] = { +/* { name, has_arg, flag, val }, where: + * name - name of long parameter + * has_arg = 0 - no argument, 1 - need argument, 2 - unnesessary argument + * flag = NULL to return val, pointer to int - to set it + * value of val (function returns 0) + * val - getopt_long return value or value, flag setting to + * !!! last string - for zeros !!! + */ + {"prefocal", 1, 0, '0'}, + {"postfocal", 1, 0, '1'}, + {"distance", 1, 0, 'D'}, + {"pixsize", 1, 0, '0'}, + { 0, 0, 0, 0 } + }; + if(argc == 1){ + usage(_("Parameters required")); + } + while(1){ + int opt; + if((opt = getopt_long(argc, argv, short_options, + long_options, NULL)) == -1) break; + switch(opt){ + /*case 0: // only long option + // do something? + break;*/ + case '0': + name0 = strdup(optarg); + break; + case '1': + name1 = strdup(optarg); + break; + case 'D': + if(!myatod(&distance, optarg)) + usage("Parameter should be a floating point number!"); + break; + case 'p': + if(!myatod(&pixsize, optarg)) + usage("Parameter should be an int or floating point number!"); + break; + default: + usage(NULL); + } + } + argc -= optind; + argv += optind; + if(argc > 0){ + // "Игнорирую аргумент[ы]:\n" + printf(_("Ignore argument[s]:\n")); + for (i = 0; i < argc; i++) + printf("%s ", argv[i]); + printf("\n"); + } + if(!name0 || !name1) + usage("You should point to both spots-files"); + if(distance < 0.) + usage("What is the distance between images?"); +} + +typedef struct{ + char *data; + size_t len; +} mmapbuf; +/** + * Mmap file to a memory area + * + * @param filename (i) - name of file to mmap + * @return stuct with mmap'ed file or die + */ +mmapbuf *My_mmap(char *filename){ + int fd; + char *ptr; + size_t Mlen; + struct stat statbuf; + if(!filename) ERRX(_("No filename given!")); + if((fd = open(filename, O_RDONLY)) < 0) + ERR(_("Can't open %s for reading"), filename); + if(fstat (fd, &statbuf) < 0) + ERR(_("Can't stat %s"), filename); + Mlen = statbuf.st_size; + if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + ERR(_("Mmap error for input")); + if(close(fd)) ERR(_("Can't close mmap'ed file")); + mmapbuf *ret = MALLOC(mmapbuf, 1); + ret->data = ptr; + ret->len = Mlen; + return ret; +} +void My_munmap(mmapbuf *b){ + if(munmap(b->data, b->len)) + ERR(_("Can't munmap")); + FREE(b); +} + +/** + * Read spots-file and fill hartmann structure + * @param filename - name of spots-file + * @return dynamically allocated hartmanogram structure + */ +hartmann *read_spots(char *filename){ + assert(filename); + mmapbuf *M = NULL; + int L = 0; + List *spots = NULL, *curspot = NULL; + M = My_mmap(filename); + hartmann *H = MALLOC(hartmann, 1); + H->filename = strdup(filename); + char *pos = M->data, *epos = pos + M->len; + for(; pos && pos < epos; pos = strchr(pos+1, '\n')){ + spot *Spot = MALLOC(spot, 1); + double x, y; + if(3 != sscanf(pos, "%d %*s %*s %*s %*s %lf %lf", &Spot->id, &x, &y)) + continue; + Spot->x = x + AxisX - ImHeight; + Spot->y = y - AxisY + ImHeight; + L++; + if(spots) + curspot = list_add(&curspot, LIST_T(Spot)); + else + curspot = list_add(&spots, LIST_T(Spot)); + }; + H->len = L; + H->spots = spots; + My_munmap(M); + return H; +} + +void h_free(hartmann **H){ + listfree_function(free); + list_free(&((*H)->spots)); + listfree_function(NULL); + free((*H)->filename); + FREE(*H); +} + +// temporary structure for building of coordinates-gradients list +typedef struct{ + double x; + double y; + double Dx; + double Dy; +} CG; +/** + * + * @param H (i) - array of thwo hartmannograms (0 - prefocal, 1 - postfocal) + * @param coords (o) - array of gradients' coordinates on prefocal image (allocated here) - the same as H[0]->spots coordinates + * @param grads (o) - gradients' array (allocated here) + * @param scale (o) - scale of polar coordinate R (== Rmax) + * @return size of built arrays + */ +size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale){ + size_t Sz = 0, i; + assert(H); assert(H[0]); assert(H[1]); + List *S0 = H[0]->spots, *S1; + List *CG_list = NULL, *curCG = NULL; + //printf(RED "\nspots\n" OLDCOLOR "\n"); + double Scale = pixsize * 1e-6 / distance / 2., S = 0.; // tg(2a)=dx/D -> a \approx dx/(2D) + /* + * Both lists have sorted structure + * but they can miss some points - that's why we should find exact points. + * To store dinamycally data I use List + */ + for(; S0; S0 = S0->next){ + spot *Sp0 = (spot*)S0->data; + int Id0 = Sp0->id; + S1 = H[1]->spots; + for(; S1; S1 = S1->next){ + spot *Sp1 = (spot*)S1->data; + if(Sp1->id > Id0) break; // point with Id0 not found + if(Sp1->id == Id0){ + CG *cg = MALLOC(CG, 1); +/* printf("id=%d (%g, %g), dx=%g, dy=%g\n", Sp0->id, Sp0->x, Sp0->y, + (Sp1->x-Sp0->x)*Scale, (Sp1->y-Sp0->y)*Scale); +*/ cg->x = Sp0->x; cg->y = Sp0->y; + cg->Dx = -(Sp1->x-Sp0->x)*Scale; + cg->Dy = -(Sp1->y-Sp0->y)*Scale; + Sz++; + if(CG_list) + curCG = list_add(&curCG, cg); + else + curCG = list_add(&CG_list, cg); + break; + } + } + } + polar *C = MALLOC(polar, Sz), *cptr = C; + point *G = MALLOC(point, Sz), *gptr = G; + curCG = CG_list; + for(i = 0; i < Sz; i++, cptr++, gptr++, curCG = curCG->next){ + double x, y, length, R; + CG *cur = (CG*)curCG->data; + x = cur->x; y = cur->y; + R = sqrt(x*x + y*y); + cptr->r = R; + if(S < R) S = R; // find max R + cptr->theta = atan2(y, x); + x = cur->Dx; y = cur->Dy; + length = sqrt(1. + x*x + y*y); // length of vector for norm + gptr->x = x / length; + gptr->y = y / length; + } + cptr = C; + for(i = 0; i < Sz; i++, cptr++) + cptr->r /= S; + *scale = S; + *coords = C; *grads = G; + listfree_function(free); + list_free(&CG_list); + listfree_function(NULL); + return Sz; +} + +int main(int argc, char **argv){ + int _U_ i; + double scale; + hartmann _U_ *images[2]; + parse_args(argc, argv); + images[0] = read_spots(name0); + images[1] = read_spots(name1); + polar *coords = NULL; point *grads = NULL; + size_t _U_ L = get_gradients(images, &coords, &grads, &scale); + h_free(&images[0]); + h_free(&images[1]); +/* printf(GREEN "\nSpots:\n" OLDCOLOR "\n\tr\ttheta\tDx\tDy\n"); + for(i = 0; i < L; i++){ + printf("%8.1f%8.4f%8.4f%8.4f\n", coords[i].r, coords[i].theta, + grads[i].x, grads[i].y); + } +*/ int Zsz, lastidx; + double *Zidxs = LS_gradZdecomposeR(15, L, coords, grads, &Zsz, &lastidx); + lastidx++; + printf("\n" RED "GradZ decompose, coefficients (%d):" OLDCOLOR "\n", lastidx); + for(i = 0; i < lastidx; i++) printf("%5.3f, ", Zidxs[i]); + printf("\n\n"); + const int GridSize = 15; + int G2 = GridSize * GridSize; + polar *rect = MALLOC(polar, G2), *rptr = rect; + double *mirZ = MALLOC(double, G2); + #define ADD 0. + int j; double Stp = 2./((double)GridSize - 1.); + for(j = 0; j < GridSize; j++){ + double y = ((double) j + ADD) * Stp - 1.; + for(i = 0; i < GridSize; i++, rptr++){ + double x = ((double) i + ADD)* Stp - 1.; + double R2 = x*x + y*y; + rptr->r = sqrt(R2); + rptr->theta = atan2(y, x); + //printf("x=%g, y=%g, r=%g, t=%g\n", x,y,rptr->r, rptr->theta); + if(R2 > 1.) continue; + mirZ[j*GridSize+i] = R2 / 32.; // mirror surface, z = r^2/(4f), or (z/R) = (r/R)^2/(4[f/R]) + //mirZ[j*GridSize+i] = R2; + } + } + printf("\n\nCoeff: %g\n\n", Zidxs[4]*sqrt(3.)); +Zidxs[4] = 0.; Zidxs[1] = 0.; Zidxs[2] = 0.; + double *comp = ZcomposeR(lastidx, Zidxs, G2, rect); + printf("\n"); + double SS = 0.; int SScntr = 0; + for(j = GridSize-1; j > -1; j--){ + for(i = 0; i < GridSize; i++){ + int idx = j*GridSize+i; + double Diff = mirZ[idx] - comp[idx]; + // printf("%7.3f", Diff); + printf("%7.3f", comp[idx]*1e3); + if(rect[idx].r < 1.){ SS += Diff; SScntr++;} + //printf("%7.3f", (comp[idx] + 0.03)/ mirZ[idx]); + } + printf("\n"); + } + SS /= SScntr; + printf("\naver: %g\n", SS); + for(j = GridSize-1; j > -1; j--){ + for(i = 0; i < GridSize; i++){ + int idx = j*GridSize+i; + //printf("%7.3f", (comp[idx] + SS) / mirZ[idx]); + double Z = (fabs(comp[idx]) < 2*DBL_EPSILON) ? 0. : comp[idx] + SS - mirZ[idx]; + printf("%7.3f", Z * 1e3); + } + printf("\n"); + } + FREE(comp); FREE(Zidxs); + FREE(coords); + FREE(grads); + return 0; +} + + diff --git a/Zernike/simple_list.c b/Zernike/simple_list.c new file mode 100644 index 0000000..cd507d3 --- /dev/null +++ b/Zernike/simple_list.c @@ -0,0 +1,112 @@ +/* + * simple_list.c - simple one-direction list + * + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include +#include + +#include "simple_list.h" + +#define MALLOC alloc_simple_list +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +static void *alloc_simple_list(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} + +void (*listdata_free)(void *node_data) = NULL; // external function to free(listdata) + +#define FREE(arg) do{free(arg); arg = NULL;}while(0) + +/** + * add element v to list with root root (also this can be last element) + * @param root (io) - pointer to root (or last element) of list or NULL + * if *root == NULL, just created node will be placed there + * @param v - data inserted + * @return pointer to created node + */ +List *list_add(List **root, void *v){ + List *node, *last; + if((node = (List*) MALLOC(1, sizeof(List))) == 0) return NULL; // allocation error + node->data = v; // insert data + if(root){ + if(*root){ // there was root node - search last + last = *root; + while(last->next) last = last->next; + last->next = node; // insert pointer to new node into last element in list + } + if(!*root) *root = node; + } + return node; +} + +/** + * set listdata_free() + * @param fn - function that freec memory of (node) + */ +void listfree_function(void (*fn)(void *node)){ + listdata_free = fn; +} + +/** + * remove all nodes in list + * @param root - pointer to root node + */ +void list_free(List **root){ + List *node = *root, *next; + if(!root || !*root) return; + do{ + next = node->next; + if(listdata_free) + listdata_free(node->data); + free(node); + node = next; + }while(node); + *root = NULL; +} + +#ifdef STANDALONE +int main(void) { + List *tp = NULL, *root_p = NULL; + int i, ins[] = {4,2,6,1,3,4,7}; + for(i = 0; i < 7; i++){ + if(!(tp = list_add(&tp, ins[i]))) + err(1, "Malloc error"); // can't insert + if(!root_p) root_p = tp; + } + tp = root_p; + i = 0; + do{ + printf("element %d = %d\n", i++, tp->data); + tp = tp->next; + }while(tp); + list_free(&root_p); + return 0; +} +#endif diff --git a/Zernike/simple_list.h b/Zernike/simple_list.h new file mode 100644 index 0000000..643ee4f --- /dev/null +++ b/Zernike/simple_list.h @@ -0,0 +1,38 @@ +/* + * simple_list.h - header file for simple list support + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __SIMPLE_LIST_H__ +#define __SIMPLE_LIST_H__ + +#define LIST_T(x) ((void*) x) + +typedef struct list_{ + void *data; + struct list_ *next; +} List; + +// add element v to list with root root (also this can be last element) +List *list_add(List **root, void *v); +// set listdata_free() +void listfree_function(void (*fn)(void *node)); +void list_free(List **root); + +#endif // __SIMPLE_LIST_H__ diff --git a/Zernike/test.c b/Zernike/test.c new file mode 100644 index 0000000..3f04e1e --- /dev/null +++ b/Zernike/test.c @@ -0,0 +1,404 @@ +/* + * test.c - test for zernike.c + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include + +#include "zernike.h" +/**/ +double IdxsOri[] = {2., // смещение, 0 + 1.1, -0.8, // наклон, 1-2 + 5.5, -3.2, 0., // астигматизм, дефокус, аст., 3-5 + 6.8, 5.5, 0., 0.,// трилистник, кома, 6-9 + 0., 0., 3.3, 1.4, 8.}; // 10-14 +/** +double IdxsOri[] = {1., // смещение, 0 + 1., 1., // наклон, 1-2 + 1., 1., 0., // астигматизм, дефокус, аст., 3-5 + 1., 1., 0., 0.,// трилистник, кома, 6-9 + 0., 0., 1., 1., 1.}; // 10-14 +**/ +const int ORI_SZ = 15; +const int W = 16, H = 16, Pmax = 8; + +//#define PLOT_ + +#define RAD 57.2957795130823 +#define D2R(x) ((x) / RAD) +#define R2D(x) ((x) * RAD) + +void plotRD(double *A, double *A1, int W, int H){ + int i,j; + double S = 0., N = 0.; + printf("\n\nRemaining abs. differences:\n"); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++){ + double d = fabs(A1[j*W+i] - A[j*W+i]); + if(d > DBL_EPSILON){ + N++; + S += d; + } + printf("%5.2f ", d); + } + printf("\n"); + } + double dA = S/N; + printf("average abs difference: %g\n", dA); + if(dA < 1.) return; + printf("Corrected diff:\n"); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++){ + double X = 0.; + if(fabs(A1[j*W+i]) > DBL_EPSILON) // X = A1[j*W+i]+dA; + X = fabs(A1[j*W+i] - A[j*W+i])-dA; + printf("%6.2f ", X); + } + printf("\n"); + } +} + +int main(int argc, char **argv){ + int i, j, lastidx; + double *Z, *Zidxs; +#if 0 + //double xstart, ystart; + double pt = 4. / ((W>H) ? W-1 : H-1); +#ifdef PLOT_ + Z = zernfunN(1, W,H, NULL); + printf("\nZernike for square matrix: \n"); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++) + printf("%5.2f ", Z[j*W+i]); + printf("\n"); + } + FREE(Z); +#endif + + Z = Zcompose(ORI_SZ, IdxsOri, W, H); +#ifdef PLOT_ + printf("\n\nDecomposition for image: \n"); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++) + printf("%6.2f ", Z[j*W+i]); + printf("\n"); + } +#endif +/* + xstart = (W-1)/2., ystart = (H-1)/2.; + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++){ + double yy = (j - ystart)*pt, xx = (i - xstart)*pt; + if((xx*xx + yy*yy) <= 1.) + Z[j*W+i] = + //100*(xx*xx*xx*xx*xx-yy*yy*yy); + 100.*(xx-0.3)*yy+200.*xx*yy*yy+300.*yy*yy*yy*yy; + printf("%6.2f ", Z[j*W+i]); + } + printf("\n"); + } +*/ + Zidxs = Zdecompose(Pmax, W, H, Z, &j, &lastidx); + printf("\nCoefficients: "); + for(i = 0; i <= lastidx; i++) printf("%5.3f, ", Zidxs[i]); + + double *Zd = Zcompose(j, Zidxs, W, H); +#ifdef PLOT_ + printf("\n\nComposed image:\n"); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++) + printf("%6.2f ", Zd[j*W+i]); + printf("\n"); + } + plotRD(Z, Zd, W, H); +#endif + FREE(Zd); + + //W--, H--; + point *testArr = MALLOC(point, W*H); + +/* + xstart = (W-1)/2., ystart = (H-1)/2.; + for(j = 0; j < H; j++){ + point *p = &testArr[j*W]; + for(i = 0; i < W; i++, p++){ + double yy = (j - ystart)*pt, xx = (i - xstart)*pt; + if((xx*xx + yy*yy) <= 1.){ + //p->x = 500.*xx*xx*xx*xx; + //p->y = -300.*yy*yy; + p->x = 100.*yy+200.*yy*yy; + p->y = 100.*(xx-0.3)+400.*xx*yy+1200.*yy*yy*yy; + } + printf("(%4.1f,%4.1f) ",p->x, p->y); + } + printf("\n"); + } +*/ + + for(j = 1; j < H-1; j++){ + point *p = &testArr[j*W+1]; + double *d = &Z[j*W+1]; + for(i = 1; i < W-1; i++, p++, d++){ + p->x = (d[1]-d[-1])/pt; + p->y = (d[W]-d[-W])/pt; + } + } +#ifdef PLOT_ + printf("\n\nGradients of field\n"); + for(j = 0; j < H; j++){ + point *p = &testArr[j*W]; + for(i = 0; i < W; i++, p++) + printf("(%4.1f,%4.1f) ",p->x, p->y); + printf("\n"); + } +#endif + //Pmax = 2; + double *_Zidxs = gradZdecompose(Pmax, W, H, testArr, &j, &lastidx); + printf("\nCoefficients: "); + for(i = 0; i <= lastidx; i++) printf("%5.3f, ", _Zidxs[i]); + //lastidx = j; + point *gZd = gradZcompose(j, _Zidxs, W, H); +#ifdef PLOT_ + printf("\n\nComposed image:\n"); + for(j = 0; j < H; j++){ + point *p = &gZd[j*W]; + for(i = 0; i < W; i++, p++) + printf("(%4.1f,%4.1f) ",p->x, p->y); + printf("\n"); + } + printf("\n\nRemaining abs differences:\n"); + for(j = 0; j < H; j++){ + point *p = &gZd[j*W]; + point *d = &testArr[j*W]; + for(i = 0; i < W; i++, p++, d++) + printf("%5.2f ",sqrt((p->x-d->x)*(p->x-d->x)+(p->y-d->y)*(p->y-d->y))); + printf("\n"); + } +#endif + FREE(gZd); + FREE(testArr); + + lastidx++; + double *ZidxsN = convGradIdxs(_Zidxs, lastidx); + printf("\nCoefficients for Zernike:\n i n m Z[i] gradZ[i] S[i] Zori[i]\n"); + for(i = 0; i < lastidx; i++){ + int n,m; + convert_Zidx(i, &n, &m); + printf("%2d%3d%3d%8.1f%8.1f%8.1f", i, n,m, Zidxs[i], ZidxsN[i],_Zidxs[i]); + if(i < ORI_SZ) printf("%8.1f", IdxsOri[i]); + printf("\n"); + } + FREE(Zidxs); +#ifdef PLOT_ + printf("\n\nComposed image:\n"); + Zd = Zcompose(lastidx, ZidxsN, W, H); + for(j = 0; j < H; j++){ + for(i = 0; i < W; i++) + printf("%6.1f ", Zd[j*W+i]); + printf("\n"); + } + plotRD(Z, Zd, W, H); + FREE(Zd); +#endif + FREE(_Zidxs); FREE(ZidxsN); +#endif // 0 + + int Sz = 256; + double dTh = D2R(360./32); + polar *P = MALLOC(polar, Sz), *Pptr = P; + void print256(double *Par){ + for(j = 0; j < 32; j++){ + for(i = 0; i < 8; i++) printf("%6.1f", Par[i*32+j]); + printf("\n"); + } + } + double Z_PREC = get_prec(); + void print256diff(double *Ori, double *App){ + printf(RED "Difference" OLDCOLOR " from original in percents:\t\t\t\t\tabs:\n"); + for(j = 0; j < 32; j++){ + for(i = 0; i < 8; i++){ + double den = Ori[i*32+j]; + if(fabs(den) > Z_PREC) + printf("%6.0f", fabs(App[i*32+j]/den-1.)*100.); + else + printf(" /0 "); + } + printf("\t"); + for(i = 0; i < 8; i++) + printf("%7.3f", fabs(App[i*32+j]-Ori[i*32+j])); + printf("\n"); + } + } + void print_std(double *p1, double *p2, int sz){ + int i; + double d = 0., d2 = 0., dmax = 0.; + for(i = 0; i < sz; i++, p1++, p2++){ + double delta = (*p2) - (*p1); + d += delta; + d2 += delta*delta; + double m = fabs(delta); + if(m > dmax) dmax = m; + } + d /= sz; + printf("\n" GREEN "Std: %g" OLDCOLOR ", max abs diff: %g\n\n", (d2/sz - d*d), dmax); + } + void print_idx_diff(double *idx){ + int i; + double *I = idx; + //printf("\nDifferences of indexes (i-i0 / i/i0):\n"); + printf("\nidx (app / real):\n"); + for(i = 0; i < ORI_SZ; i++, I++) + printf("%d: (%3.1f / %3.1f), ", i, *I, IdxsOri[i]); + //printf("%d: (%3.1f / %3.1f), ", i, *I - IdxsOri[i], *I/IdxsOri[i]); + print_std(idx, IdxsOri, ORI_SZ); + printf("\n"); + } + void mktest__(double* (*decomp_fn)(int, int, polar*, double*, int*, int*), char *msg){ + Zidxs = decomp_fn(Pmax, Sz, P, Z, &j, &lastidx); + printf("\n" RED "%s, coefficients (%d):" OLDCOLOR "\n", msg, lastidx); + lastidx++; + for(i = 0; i < lastidx; i++) printf("%5.3f, ", Zidxs[i]); + printf("\nComposing: %s(%d, Zidxs, %d, P)\n", msg, lastidx, Sz); + double *comp = ZcomposeR(lastidx, Zidxs, Sz, P); + print_std(Z, comp, Sz); + print_idx_diff(Zidxs); + print256diff(Z, comp); + FREE(comp); + FREE(Zidxs); + printf("\n"); + } + #define mktest(a) mktest__(a, #a) + double R_holes[] = {.175, .247, .295, .340, .379, .414, .448, .478}; + for(i = 0; i < 8; i++){ + // double RR = (double)i * 0.1 + 0.3; + // double RR = (double)i * 0.14+0.001; + // double RR = (double)i * 0.07 + 0.5; + double RR = R_holes[i]; // / 0.5; + double Th = 0.; + for(j = 0; j < 32; j++, Pptr++, Th += dTh){ + Pptr->r = RR; + Pptr->theta = Th; + } + } + Z = ZcomposeR(ORI_SZ, IdxsOri, Sz, P); + printf(RED "Original:" OLDCOLOR "\n"); + print256(Z); + + mktest(ZdecomposeR); + mktest(QR_decompose); + mktest(LS_decompose); + + +// ann_Z(4, Sz, P, NULL); + //polar P1[] = {{0.2, 0.}, {0.6, M_PI_2}, {0.1, M_PI}, {0.92, 3.*M_PI_2}}; + //ann_Z(8, 4, P1, NULL); +/* + double id_ann[] = {1.,2.,0.1,-1.,0.5}; + Z = ann_Zcompose(5, id_ann, Sz, P); +*/ + FREE(Z); + Z = ann_Zcompose(ORI_SZ, IdxsOri, Sz, P); + printf(RED "Annular:" OLDCOLOR "\n"); + print256(Z); + Zidxs = ann_Zdecompose(7, Sz, P, Z, &j, &lastidx); + printf("\nann_Zdecompose, coefficients:\n"); + for(i = 0; i <= lastidx; i++) printf("%5.3f, ", Zidxs[i]); + double *comp = ann_Zcompose(lastidx, Zidxs, Sz, P); + print_std(Z, comp, Sz); + print_idx_diff(Zidxs); + print256diff(Z, comp); + FREE(comp); + FREE(Zidxs); + +/* + * TEST for gradients on hartmann's net + */ + + point *gradZ = MALLOC(point, Sz); + Pptr = P; + point *dP = gradZ; + for(i = 0; i < Sz; i++, Pptr++, dP++){ + double x = Pptr->r * cos(Pptr->theta), y = Pptr->r * sin(Pptr->theta); + double x2 _U_ = x*x, x3 _U_ = x*x2; + double y2 _U_ = y*y, y3 _U_ = y*y2; + double sx, cx; + sincos(x, &sx, &cx); + + double val = 1000.*x3 + y2 + 30.*cx*y3; + Z[i] = val; val *= 0.05; + dP->x = 3000.*x2 - 30.*sx*y3 + val * drand48(); + dP->y = 2.*y + 90.*cx*y2 + val * drand48(); + /* + double val = 50.*x3*y3 + 3.*x2*y2 - 2.*x*y3 - 8.*x3*y + 7.*x*y; + Z[i] = val; val *= 0.05; // 5% from value + dP->x = 150.*x2*y3 + 6.*x*y2 - 2.*y3 - 24.*x2*y + 7.*y + val * drand48(); // + error 5% + dP->y = 150.*x3*y2 + 6.*x2*y - 6.*x*y2 - 8.*x3 + 7.*x + val * drand48(); + */ + } + printf("\n" RED "Z cubic:" OLDCOLOR "\n"); + print256(Z); + mktest(ZdecomposeR); + mktest(LS_decompose); + mktest(QR_decompose); + + Zidxs = LS_gradZdecomposeR(Pmax+2, Sz, P, gradZ, &i, &lastidx); + lastidx++; + printf("\n" RED "GradZ decompose, coefficients (%d):" OLDCOLOR "\n", lastidx); + for(i = 0; i < lastidx; i++) printf("%5.3f, ", Zidxs[i]); + comp = ZcomposeR(lastidx, Zidxs, Sz, P); + double averD = 0.; + for(i = 0; i < Sz; i++) averD += Z[i] - comp[i]; + averD /= (double)Sz; + printf("Z0 = %g\n", averD); + for(i = 0; i < Sz; i++) comp[i] += averD; + print256diff(Z, comp); + FREE(comp); + FREE(Zidxs); + + Zidxs = gradZdecomposeR(Pmax+2, Sz, P, gradZ, &i, &lastidx); +/* + point *gZd = gradZcomposeR(lastidx, Zidxs, Sz, P); + printf("\n\nRemaining abs differences:\n"); + for(i = 0; i < Sz; i++){ + point p = gZd[i]; + point d = gradZ[i]; + printf("%5.2f ",sqrt((p.x-d.x)*(p.x-d.x)+(p.y-d.y)*(p.y-d.y))); + } + printf("\n"); + FREE(gZd); +*/ +// double *ZidxsN = convGradIdxs(Zidxs, lastidx); +// FREE(Zidxs); Zidxs = ZidxsN; + lastidx++; + printf("\n" RED "GradZ decompose, coefficients (%d):" OLDCOLOR "\n", lastidx); + for(i = 0; i < lastidx; i++) printf("%5.3f, ", Zidxs[i]); + comp = ZcomposeR(lastidx, Zidxs, Sz, P); + averD = 0.; + for(i = 0; i < Sz; i++) averD += Z[i] - comp[i]; + averD /= (double)Sz; + printf("Z0 = %g\n", averD); + for(i = 0; i < Sz; i++) comp[i] += averD; + print256diff(Z, comp); + FREE(comp); + FREE(Zidxs); + + FREE(Z); + return 0; +} diff --git a/Zernike/zern_private.h b/Zernike/zern_private.h new file mode 100644 index 0000000..dd710be --- /dev/null +++ b/Zernike/zern_private.h @@ -0,0 +1,52 @@ +/* + * zern_private.h - private variables for zernike.c & zernikeR.c + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __ZERN_PRIVATE_H__ +#define __ZERN_PRIVATE_H__ + +#include +#include +#include +#include + +#ifndef iabs +#define iabs(a) (((a)<(0)) ? (-a) : (a)) +#endif + +#ifndef DBG +#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0) +#endif + +extern double *FK; +extern double Z_prec; + +// zernike.c +void build_factorial(); +void free_rpow(double ***Rpow, int n); +void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow); +double **build_rpowR(int n, int Sz, polar *P); + + + +// zernike_annular.c +polar *conv_r(polar *r0, int Sz); + +#endif // __ZERN_PRIVATE_H__ diff --git a/Zernike/zernike.c b/Zernike/zernike.c new file mode 100644 index 0000000..1fcf872 --- /dev/null +++ b/Zernike/zernike.c @@ -0,0 +1,515 @@ +/* + * zernike.c - wavefront decomposition using Zernike polynomials + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +/* + * Literature: +@ARTICLE{2007OExpr..1518014Z, + author = {{Zhao}, C. and {Burge}, J.~H.}, + title = "{Orthonormal vector polynomials in a unit circle Part I: basis set derived from gradients of Zernike polynomials}", + journal = {Optics Express}, + year = 2007, + volume = 15, + pages = {18014}, + doi = {10.1364/OE.15.018014}, + adsurl = {http://adsabs.harvard.edu/abs/2007OExpr..1518014Z}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} +@ARTICLE{2008OExpr..16.6586Z, + author = {{Zhao}, C. and {Burge}, J.~H.}, + title = "{Orthonormal vector polynomials in a unit circle, Part II : completing the basis set}", + journal = {Optics Express}, + year = 2008, + month = apr, + volume = 16, + pages = {6586}, + doi = {10.1364/OE.16.006586}, + adsurl = {http://adsabs.harvard.edu/abs/2008OExpr..16.6586Z}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + * + * !!!ATTENTION!!! Axe Y points upside-down! + */ + + +#include "zernike.h" +#include "zern_private.h" + +/** + * safe memory allocation for macro ALLOC + * @param N - number of elements to allocate + * @param S - size of single element (typically sizeof) + * @return pointer to allocated memory area + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} + +double *FK = NULL; +/** + * Build pre-computed array of factorials from 1 to 100 + */ +void build_factorial(){ + double F = 1.; + int i; + if(FK) return; + FK = MALLOC(double, 100); + FK[0] = 1.; + for(i = 1; i < 100; i++) + FK[i] = (F *= (double)i); +} + +double Z_prec = 1e-6; // precision of Zernike coefficients +/** + * Set precision of Zernice coefficients decomposition + * @param val - new precision value + */ +void set_prec(double val){ + Z_prec = val; +} + +double get_prec(){ + return Z_prec; +} + +/** + * Convert polynomial order in Noll notation into n/m + * @param p (i) - order of Zernike polynomial in Noll notation + * @param N (o) - order of polynomial + * @param M (o) - angular parameter + */ +void convert_Zidx(int p, int *N, int *M){ + int n = (int) floor((-1.+sqrt(1.+8.*p)) / 2.); + *M = (int)(2.0*(p - n*(n+1.)/2. - 0.5*(double)n)); + *N = n; +} + +/** + * Free array of R powers with power n + * @param Rpow (i) - array to free + * @param n - power of Zernike polinomial for that array (array size = n+1) + */ +void free_rpow(double ***Rpow, int n){ + int i, N = n+1; + for(i = 0; i < N; i++) FREE((*Rpow)[i]); + FREE(*Rpow); +} + +/** + * Build two arrays: with R and its powers (from 0 to n inclusive) + * for cartesian quoter I of matrix with size WxH + * @param W, H - size of initial matrix + * @param n - power of Zernike polinomial (array size = n+1) + * @param Rad (o) - NULL or array with R in quater I + * @param Rad_pow (o) - NULL or array with powers of R + */ +void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow){ + double Rnorm = fmax((W - 1.) / 2., (H - 1.) / 2.); + int i,j, k, N = n+1, w = (W+1)/2, h = (H+1)/2, S = w*h; + double **Rpow = MALLOC(double*, N); // powers of R + Rpow[0] = MALLOC(double, S); + for(j = 0; j < S; j++) Rpow[0][j] = 1.; // zero's power + double *R = MALLOC(double, S); + // first - fill array of R + double xstart = (W%2) ? 0. : 0.5, ystart = (H%2) ? 0. : 0.5; + for(j = 0; j < h; j++){ + double *pt = &R[j*w], x, ww = w; + for(x = xstart; x < ww; x++, pt++){ + double yy = (j + ystart)/Rnorm, xx = x/Rnorm; + *pt = sqrt(xx*xx+yy*yy); + } + } + for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical) + Rpow[i] = MALLOC(double, S); + double *rp = Rpow[i], *rpo = Rpow[i-1]; + for(j = 0; j < h; j++){ + int idx = j*w; + double *pt = &rp[idx], *pto = &rpo[idx], *r = &R[idx]; + for(k = 0; k < w; k++, pt++, pto++, r++) + *pt = (*pto) * (*r); // R^{n+1} = R^n * R + } + } + if(Rad) *Rad = R; + else FREE(R); + if(Rad_pow) *Rad_pow = Rpow; + else free_rpow(&Rpow, n); +} + +/** + * Calculate value of Zernike polynomial on rectangular matrix WxH pixels + * Center of matrix will be zero point + * Scale will be set by max(W/2,H/2) + * @param n - order of polynomial (max: 100!) + * @param m - angular parameter of polynomial + * @param W - width of output array + * @param H - height of output array + * @param norm (o) - (if !NULL) normalize factor + * @return array of Zernike polynomials on given matrix + */ +double *zernfun(int n, int m, int W, int H, double *norm){ + double Z = 0., *Zarr = NULL; + bool erparm = false; + if(W < 2 || H < 2) + errx(1, "Sizes of matrix must be > 2!"); + if(n > 100) + errx(1, "Order of Zernike polynomial must be <= 100!"); + if(n < 0) erparm = true; + if(n < iabs(m)) erparm = true; // |m| must be <= n + int d = n - m; + if(d % 2) erparm = true; // n-m must differ by a prod of 2 + if(erparm) + errx(1, "Wrong parameters of Zernike polynomial (%d, %d)", n, m); + if(!FK) build_factorial(); + double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2.; // coordinate of circle's middle + int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, w = (W+1)/2; + double *R, **Rpow; + build_rpow(W, H, n, &R, &Rpow); + int SS = W * H; + double ZSum = 0.; + // now fill output matrix + Zarr = MALLOC(double, SS); // output matrix W*H pixels + for(j = 0; j < H; j++){ + double *Zptr = &Zarr[j*W]; + double Ryd = fabs(j - Yc); + int Ry = w * (int)Ryd; // Y coordinate on R matrix + for(i = 0; i < W; i++, Zptr++){ + Z = 0.; + double Rxd = fabs(i - Xc); + int Ridx = Ry + (int)Rxd; // coordinate on R matrix + if(R[Ridx] > 1.) continue; // throw out points with R>1 + // calculate R_n^m + for(k = 0; k <= iup; k++){ // Sum + double z = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[(n+m_abs)/2-k]*FK[(n-m_abs)/2-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); + Z += z * Rpow[n-2*k][Ridx]; // *R^{n-2k} + } + // normalize + double eps_m = (m) ? 1. : 2.; + Z *= sqrt(2.*(n+1.) / M_PI / eps_m); + double theta = atan2(j - Yc, i - Xc); + // multiply to angular function: + if(m){ + if(m > 0) + Z *= cos(theta*(double)m_abs); + else + Z *= sin(theta*(double)m_abs); + } + *Zptr = Z; + ZSum += Z*Z; + } + } + if(norm) *norm = ZSum; + // free unneeded memory + FREE(R); + free_rpow(&Rpow, n); + return Zarr; +} + +/** + * Zernike polynomials in Noll notation for square matrix + * @param p - index of polynomial + * @other params are like in zernfun + * @return zernfun + */ +double *zernfunN(int p, int W, int H, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + return zernfun(n,m,W,H,norm); +} + +/** + * Zernike decomposition of image 'image' WxH pixels + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param W, H (i) - size of image + * @param image(i) - image itself + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx (o) - (if !NULL) last non-zero coefficient + * @return array of Zernike coefficients + */ +double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx){ + int i, SS = W*H, pmax, maxIdx = 0; + double *Zidxs = NULL, *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(double, W*H); + memcpy(icopy, image, W*H*sizeof(double)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 0; i < pmax; i++){ // now we fill array + double norm, *Zcoeff = zernfunN(i, W, H, &norm); + int j; + double *iptr = icopy, *zptr = Zcoeff, K = 0.; + for(j = 0; j < SS; j++, iptr++, zptr++) + K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient + Zidxs[i] = K; + if(fabs(K) < Z_prec){ + Zidxs[i] = 0.; + continue; // there's no need to substract values that are less than our precision + } + maxIdx = i; + iptr = icopy; zptr = Zcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++) + *iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values + FREE(Zcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Zernike restoration of image WxH pixels by Zernike polynomials coefficients + * @param Zsz (i) - number of elements in coefficients array + * @param Zidxs(i) - array with Zernike coefficients + * @param W, H (i) - size of image + * @return restored image + */ +double *Zcompose(int Zsz, double *Zidxs, int W, int H){ + int i, SS = W*H; + double *image = MALLOC(double, SS); + for(i = 0; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + double *Zcoeff = zernfunN(i, W, H, NULL); + int j; + double *iptr = image, *zptr = Zcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++) + *iptr += K * (*zptr); // add next Zernike polynomial + FREE(Zcoeff); + } + return image; +} + + +/** + * Components of Zj gradient without constant factor + * all parameters are as in zernfun + * @return array of gradient's components + */ +point *gradZ(int n, int m, int W, int H, double *norm){ + point *gZ = NULL; + bool erparm = false; + if(W < 2 || H < 2) + errx(1, "Sizes of matrix must be > 2!"); + if(n > 100) + errx(1, "Order of gradient of Zernike polynomial must be <= 100!"); + if(n < 1) erparm = true; + if(n < iabs(m)) erparm = true; // |m| must be <= n + int d = n - m; + if(d % 2) erparm = true; // n-m must differ by a prod of 2 + if(erparm) + errx(1, "Wrong parameters of gradient of Zernike polynomial (%d, %d)", n, m); + if(!FK) build_factorial(); + double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2.; // coordinate of circle's middle + int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2, w = (W+1)/2; + double *R, **Rpow; + build_rpow(W, H, n, &R, &Rpow); + int SS = W * H; + // now fill output matrix + gZ = MALLOC(point, SS); + double ZSum = 0.; + for(j = 0; j < H; j++){ + point *Zptr = &gZ[j*W]; + double Ryd = fabs(j - Yc); + int Ry = w * (int)Ryd; // Y coordinate on R matrix + for(i = 0; i < W; i++, Zptr++){ + double Rxd = fabs(i - Xc); + int Ridx = Ry + (int)Rxd; // coordinate on R matrix + double Rcurr = R[Ridx]; + if(Rcurr > 1. || fabs(Rcurr) < DBL_EPSILON) continue; // throw out points with R>1 + double theta = atan2(j - Yc, i - Xc); + // components of grad Zj: + + // 1. Theta_j + double Tj = 1., costm, sintm; + sincos(theta*(double)m_abs, &sintm, &costm); + if(m){ + if(m > 0) + Tj = costm; + else + Tj = sintm; + } + + // 2. dTheta_j/Dtheta + double dTj = 0.; + if(m){ + if(m < 0) + dTj = m_abs * costm; + else + dTj = -m_abs * sintm; + } + // 3. R_j & dR_j/dr + double Rj = 0., dRj = 0.; + for(k = 0; k <= iup; k++){ + double rj = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[isum-k]*FK[iup-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); +//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k); + int kidx = n - 2*k; + Rj += rj * Rpow[kidx][Ridx]; // *R^{n-2k} + if(kidx > 0) + dRj += rj * kidx * Rpow[kidx - 1][Ridx]; + } + // normalisation for Zernike + double eps_m = (m) ? 1. : 2., sq = sqrt(2.*(double)(n+1) / M_PI / eps_m); + Rj *= sq, dRj *= sq; + // 4. sin/cos + double sint, cost; + sincos(theta, &sint, &cost); + + // projections of gradZj + double TdR = Tj * dRj, RdT = Rj * dTj / Rcurr; + Zptr->x = TdR * cost - RdT * sint; + Zptr->y = TdR * sint + RdT * cost; + // norm factor + ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y; + } + } + if(norm) *norm = ZSum; + // free unneeded memory + FREE(R); + free_rpow(&Rpow, n); + return gZ; +} + +/** + * Build components of vector polynomial Sj + * all parameters are as in zernfunN + * @return Sj(n,m) on image points + */ +point *zerngrad(int p, int W, int H, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + if(n < 1) errx(1, "Order of gradient Z must be > 0!"); + int m_abs = iabs(m); + int i,j; + double K = 1., K1 = 1.;///sqrt(2.*n*(n+1.)); + point *Sj = MALLOC(point, W*H); + point *Zj = gradZ(n, m, W, H, NULL); + double Zsum = 0.; + if(n == m_abs || n < 3){ // trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0 + for(j = 0; j < H; j++){ + int idx = j*W; + point *Sptr = &Sj[idx], *Zptr = &Zj[idx]; + for(i = 0; i < W; i++, Sptr++, Zptr++){ + Sptr->x = K1*Zptr->x; + Sptr->y = K1*Zptr->y; + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + } + }else{ + K = sqrt(((double)n+1.)/(n-1.)); + //K1 /= sqrt(2.); + // n != |m| + // I run gradZ() twice! But another variant (to make array of Zj) can meet memory lack + point *Zj_= gradZ(n-2, m, W, H, NULL); + for(j = 0; j < H; j++){ + int idx = j*W; + point *Sptr = &Sj[idx], *Zptr = &Zj[idx], *Z_ptr = &Zj_[idx]; + for(i = 0; i < W; i++, Sptr++, Zptr++, Z_ptr++){ + Sptr->x = K1*(Zptr->x - K * Z_ptr->x); + Sptr->y = K1*(Zptr->y - K * Z_ptr->y); + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + } + FREE(Zj_); + } + FREE(Zj); + if(norm) *norm = Zsum; + return Sj; +} + +/** + * Decomposition of image with normals to wavefront by Zhao's polynomials + * all like Zdecompose + * @return array of coefficients + */ +double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx){ + int i, SS = W*H, pmax, maxIdx = 0; + double *Zidxs = NULL; + point *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(point, SS); + memcpy(icopy, image, SS*sizeof(point)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 1; i < pmax; i++){ // now we fill array + double norm; + point *dZcoeff = zerngrad(i, W, H, &norm); + int j; + point *iptr = icopy, *zptr = dZcoeff; + double K = 0.; + for(j = 0; j < SS; j++, iptr++, zptr++) + K += (zptr->x*iptr->x + zptr->y*iptr->y) / norm; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = dZcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++){ + iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values + iptr->y -= K * zptr->y; + } + FREE(dZcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of image with normals Zhao's polynomials coefficients + * all like Zcompose + * @return restored image + */ +point *gradZcompose(int Zsz, double *Zidxs, int W, int H){ + int i, SS = W*H; + point *image = MALLOC(point, SS); + for(i = 1; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + point *Zcoeff = zerngrad(i, W, H, NULL); + int j; + point *iptr = image, *zptr = Zcoeff; + for(j = 0; j < SS; j++, iptr++, zptr++){ + iptr->x += K * zptr->x; + iptr->y += K * zptr->y; + } + FREE(Zcoeff); + } + return image; +} + +double *convGradIdxs(double *gradIdxs, int Zsz){ + double *idxNew = MALLOC(double, Zsz); + int i; + for(i = 1; i < Zsz; i++){ + int n,m; + convert_Zidx(i, &n, &m); + int j = ((n+2)*(n+4) + m) / 2; + if(j >= Zsz) j = 0; + idxNew[i] = (gradIdxs[i] - sqrt((n+3.)/(n+1.))*gradIdxs[j]); + } + return idxNew; +} diff --git a/Zernike/zernike.h b/Zernike/zernike.h new file mode 100644 index 0000000..6bc9feb --- /dev/null +++ b/Zernike/zernike.h @@ -0,0 +1,92 @@ +/* + * zernike.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __ZERNIKE_H__ +#define __ZERNIKE_H__ + +#include +#include + +#define _U_ __attribute__((unused)) +#define RED "\033[1;31;40m" +#define GREEN "\033[1;32;40m" +#define OLDCOLOR "\033[0;0;0m" + +/*************** Data structures & typedefs ***************/ +// point coordinates +typedef struct{ + double x,y; +} point; + +typedef struct{ + double r,theta; +} polar; + +// 2D array +typedef struct{ + double **data; + size_t len; // size of 1D arrays + size_t num; // number of 1D arrays +}_2D; + +#ifndef DBL_EPSILON +#define DBL_EPSILON 2.2204460492503131e-16 +#endif + +/*************** Memory allocation ***************/ +#define MALLOC(type, size) ((type*)my_alloc(size, sizeof(type))) +#define FREE(ptr) do{free(ptr); ptr = NULL;}while(0) +void *my_alloc(size_t N, size_t S); + +/*************** Base functions ***************/ +void convert_Zidx(int p, int *N, int *M); +void set_prec(double val); +double get_prec(); + +/*************** Zernike on rectangular equidistant coordinate matrix ***************/ +double *zernfun(int n, int m, int W, int H, double *norm); +double *zernfunN(int p, int W, int H, double *norm); + +double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx); +double *Zcompose(int Zsz, double *Zidxs, int W, int H); + +double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx); +point *gradZcompose(int Zsz, double *Zidxs, int W, int H); +double *convGradIdxs(double *gradIdxs, int Zsz); + +/*************** Zernike on a points set ***************/ +double *zernfunR(int n, int m, int Sz, polar *P, double *norm); +double *zernfunNR(int p, int Sz, polar *P, double *norm); +double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P); +double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); +double *gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx); +double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx); +point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P); + +/*************** Annular Zernike ***************/ +_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm); +double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P); +double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx); + +#endif // __ZERNIKE_H__ diff --git a/Zernike/zernikeR.c b/Zernike/zernikeR.c new file mode 100644 index 0000000..7f43dd1 --- /dev/null +++ b/Zernike/zernikeR.c @@ -0,0 +1,616 @@ +/* + * zernikeR.c - Zernike decomposition for scattered points & annular aperture + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include + +#include "zernike.h" +#include "zern_private.h" + +/** + * Build array with R powers (from 0 to n inclusive) + * @param n - power of Zernike polinomial (array size = n+1) + * @param Sz - size of P array + * @param P (i) - polar coordinates of points + */ +double **build_rpowR(int n, int Sz, polar *P){ + int i, j, N = n + 1; + double **Rpow = MALLOC(double*, N); + Rpow[0] = MALLOC(double, Sz); + for(i = 0; i < Sz; i++) Rpow[0][i] = 1.; // zero's power + for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical) + Rpow[i] = MALLOC(double, Sz); + double *rp = Rpow[i], *rpo = Rpow[i-1]; + polar *p = P; + for(j = 0; j < Sz; j++, rp++, rpo++, p++){ + *rp = (*rpo) * p->r; + } + } + return Rpow; +} + +bool check_parameters(n, m, Sz, P){ + bool erparm = false; + if(Sz < 3 || !P) + errx(1, "Size of matrix must be > 2!"); + if(n > 100) + errx(1, "Order of Zernike polynomial must be <= 100!"); + if(n < 0) erparm = true; + if(n < iabs(m)) erparm = true; // |m| must be <= n + if((n - m) % 2) erparm = true; // n-m must differ by a prod of 2 + if(erparm) + errx(1, "Wrong parameters of Zernike polynomial (%d, %d)", n, m); + else + if(!FK) build_factorial(); + return erparm; +} + +/** + * Zernike function for scattering data + * @param n,m - orders of polynomial + * @param Sz - number of points + * @param P(i) - array with points coordinates (polar, r<=1) + * @param norm(o) - (optional) norm coefficient + * @return dynamically allocated array with Z(n,m) for given array P + */ +double *zernfunR(int n, int m, int Sz, polar *P, double *norm){ + if(check_parameters(n, m, Sz, P)) return NULL; + int j, k, m_abs = iabs(m), iup = (n-m_abs)/2; + double **Rpow = build_rpowR(n, Sz, P); + double ZSum = 0.; + // now fill output matrix + double *Zarr = MALLOC(double, Sz); // output matrix + double *Zptr = Zarr; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + double Z = 0.; + if(p->r > 1.) continue; // throw out points with R>1 + // calculate R_n^m + for(k = 0; k <= iup; k++){ // Sum + double z = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[(n+m_abs)/2-k]*FK[(n-m_abs)/2-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); + Z += z * Rpow[n-2*k][j]; // *R^{n-2k} + } + // normalize + double eps_m = (m) ? 1. : 2.; + Z *= sqrt(2.*(n+1.) / M_PI / eps_m ); + double m_theta = (double)m_abs * p->theta; + // multiply to angular function: + if(m){ + if(m > 0) + Z *= cos(m_theta); + else + Z *= sin(m_theta); + } + *Zptr = Z; + ZSum += Z*Z; + } + if(norm) *norm = ZSum; + // free unneeded memory + free_rpow(&Rpow, n); + return Zarr; +} + +/** + * Zernike polynomials in Noll notation for square matrix + * @param p - index of polynomial + * @other params are like in zernfunR + * @return zernfunR + */ +double *zernfunNR(int p, int Sz, polar *P, double *norm){ + int n, m; + convert_Zidx(p, &n, &m); + return zernfunR(n,m,Sz,P,norm); +} + +/** + * Zernike decomposition of image 'image' WxH pixels + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param heights(i) - wavefront walues in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of Zernike coefficients + */ +double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + double *Zidxs = NULL, *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(double, Sz); + memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 0; i < pmax; i++){ // now we fill array + double norm, *Zcoeff = zernfunNR(i, Sz, P, &norm); + int j; + double *iptr = icopy, *zptr = Zcoeff, K = 0.; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + + iptr = icopy; zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values + + FREE(Zcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of image in points P by Zernike polynomials coefficients + * @param Zsz (i) - number of actual elements in coefficients array + * @param Zidxs(i) - array with Zernike coefficients + * @param Sz, P(i) - number (Sz) of points (P) + * @return restored image + */ +double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + double *image = MALLOC(double, Sz); + for(i = 0; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + double *Zcoeff = zernfunNR(i, Sz, P, NULL); + int j; + double *iptr = image, *zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr += K * (*zptr); // add next Zernike polynomial + FREE(Zcoeff); + } + return image; +} + +/** + * Prints out GSL matrix + * @param M (i) - matrix to print + */ +void print_matrix(gsl_matrix *M){ + int x,y, H = M->size1, W = M->size2; + size_t T = M->tda; + printf("\nMatrix %dx%d:\n", W, H); + for(y = 0; y < H; y++){ + double *ptr = &(M->data[y*T]); + printf("str %6d", y); + for(x = 0; x < W; x++, ptr++) + printf("%6.1f ", *ptr); + printf("\n"); + } + printf("\n\n"); +} + +/* +To try less-squares fit I decide after reading of +@ARTICLE{2010ApOpt..49.6489M, + author = {{Mahajan}, V.~N. and {Aftab}, M.}, + title = "{Systematic comparison of the use of annular and Zernike circle polynomials for annular wavefronts}", + journal = {\ao}, + year = 2010, + month = nov, + volume = 49, + pages = {6489}, + doi = {10.1364/AO.49.006489}, + adsurl = {http://adsabs.harvard.edu/abs/2010ApOpt..49.6489M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} +*/ +/* + * n'th column of array m is polynomial of n'th degree, + * m'th row is m'th data point + * + * We fill matrix with polynomials by known datapoints coordinates, + * after that by less-square fit we get coefficients of decomposition + */ +double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation +/* + (from GSL) + typedef struct { + size_t size1; // rows (height) + size_t size2; // columns (width) + size_t tda; // data width (aligned) - data stride + double * data; // pointer to block->data (matrix data itself) + gsl_block * block; // block with matrix data (block->size is size) + int owner; // ==1 if matrix owns 'block' (then block will be freed with matrix) + } gsl_matrix; +*/ + // Now allocate matrix: its Nth row is equation for Nth data point, + // Mth column is Z_M + gsl_matrix *M = gsl_matrix_calloc(Sz, pmax); + // fill matrix with coefficients + int x,y; + size_t T = M->tda; + for(x = 0; x < pmax; x++){ + double norm, *Zcoeff = zernfunNR(x, Sz, P, &norm), *Zptr = Zcoeff; + double *ptr = &(M->data[x]); + for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial + *ptr = (*Zptr); + FREE(Zcoeff); + } + + gsl_vector *ans = gsl_vector_calloc(pmax); + + gsl_vector_view b = gsl_vector_view_array(heights, Sz); + gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M)) + gsl_linalg_QR_decomp(M, tau); + + gsl_vector *residual = gsl_vector_calloc(Sz); + gsl_linalg_QR_lssolve(M, tau, &b.vector, ans, residual); + double *ptr = ans->data; + int maxIdx = 0; + double *Zidxs = MALLOC(double, pmax); + for(x = 0; x < pmax; x++, ptr++){ + if(fabs(*ptr) < Z_prec) continue; + Zidxs[x] = *ptr; + maxIdx = x; + } + + gsl_matrix_free(M); + gsl_vector_free(ans); + gsl_vector_free(tau); + gsl_vector_free(residual); + + if(lastIdx) *lastIdx = maxIdx; + return Zidxs; +} + +/* + * Pseudo-annular Zernike polynomials + * They are just a linear composition of Zernike, made by Gram-Schmidt ortogonalisation + * Real Zernike koefficients restored by reverce transformation matrix + * + * Suppose we have a wavefront data in set of scattered points ${(x,y)}$, we want to + * find Zernike coefficients $z$ such that product of Zernike polynomials, $Z$, and + * $z$ give us that wavefront data with very little error (depending on $Z$ degree). + * + * Of cource, $Z$ isn't orthonormal on our basis, so we need to create some intermediate + * polynomials, $U$, which will be linear dependent from $Z$ (and this dependency + * should be strict and reversable, otherwise we wouldn't be able to reconstruct $z$ + * from $u$): $U = Zk$. So, we have: $W = Uu + \epsilon$ and $W = Zz + \epsilon$. + * + * $U$ is orthonormal, so $U^T = U^{-1}$ and (unlike to $Z$) this reverce matrix + * exists and is unique. This mean that $u = W^T U = U^T W$. + * Coefficients matrix, $k$ must be inversable, so $Uk^{-1} = Z$, this mean that + * $z = uk^{-1}$. + * + * Our main goal is to find that matrix $k$. + * + * 1. QR-decomposition + * In this case non-orthogonal matrix $Z$ decomposed to orthogonal matrix $Q$ and + * right-triangular matrix $R$. In our case $Q$ is $U$ itself and $R$ is $k^{-1}$. + * QR-decomposition gives us an easy way to compute coefficient's matrix, $k$. + * Polynomials in $Q$ are linear-dependent from $Z$, each $n^{th}$ polynomial in $Q$ + * depends from first $n$ polynomials in $Z$. So, columns of $R$ are coefficients + * for making $U$ from $Z$; rows in $R$ are coefficients for making $z$ from $u$. + * + * 2. Cholesky decomposition + * In this case for any non-orthogonal matrix with real values we have: + * $A^{T}A = LL^{T}$, where $L$ is left-triangular matrix. + * Orthogonal basis made by equation $U = A(L^{-1})^T$. And, as $A = UL^T$, we + * can reconstruct coefficients matrix: $k = (L^{-1})^T$. + */ + +double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + if(Sz < pmax) errx(1, "Number of points must be >= number of polynomials!"); + int k, x,y; + + //make new polar + polar *Pn = conv_r(P, Sz); + // Now allocate matrix: its Nth row is equation for Nth data point, + // Mth column is Z_M + gsl_matrix *M = gsl_matrix_calloc(Sz, pmax); + // Q-matrix (its first pmax columns is our basis) + gsl_matrix *Q = gsl_matrix_calloc(Sz, Sz); + // R-matrix (coefficients) + gsl_matrix *R = gsl_matrix_calloc(Sz, pmax); + // fill matrix with coefficients + size_t T = M->tda; + for(x = 0; x < pmax; x++){ + double norm, *Zcoeff = zernfunNR(x, Sz, Pn, &norm), *Zptr = Zcoeff; + double *ptr = &(M->data[x]); + for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial + *ptr = (*Zptr) / norm; + FREE(Zcoeff); + } + + gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M)) + gsl_linalg_QR_decomp(M, tau); + + gsl_linalg_QR_unpack(M, tau, Q, R); +//print_matrix(R); + gsl_matrix_free(M); + gsl_vector_free(tau); + + double *Zidxs = MALLOC(double, pmax); + T = Q->tda; + size_t Tr = R->tda; + + for(k = 0; k < pmax; k++){ // cycle by powers + double sumD = 0.; + double *Qptr = &(Q->data[k]); + for(y = 0; y < Sz; y++, Qptr+=T){ // cycle by points + sumD += heights[y] * (*Qptr); + } + Zidxs[k] = sumD; + } + gsl_matrix_free(Q); + + // now restore Zernike coefficients + double *Zidxs_corr = MALLOC(double, pmax); + int maxIdx = 0; + for(k = 0; k < pmax; k++){ + double sum = 0.; + double *Rptr = &(R->data[k*(Tr+1)]), *Zptr = &(Zidxs[k]); + for(x = k; x < pmax; x++, Zptr++, Rptr++){ + sum += (*Zptr) * (*Rptr); + } + if(fabs(sum) < Z_prec) continue; + Zidxs_corr[k] = sum; + maxIdx = k; + } + gsl_matrix_free(R); + FREE(Zidxs); + FREE(Pn); + if(lastIdx) *lastIdx = maxIdx; + return Zidxs_corr; +} + + +/** + * Components of Zj gradient without constant factor + * @param n,m - orders of polynomial + * @param Sz - number of points + * @param P (i) - coordinates of reference points + * @param norm (o) - norm factor or NULL + * @return array of gradient's components + */ +point *gradZR(int n, int m, int Sz, polar *P, double *norm){ + if(check_parameters(n, m, Sz, P)) return NULL; + point *gZ = NULL; + int j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2; + double **Rpow = build_rpowR(n, Sz, P); + // now fill output matrix + gZ = MALLOC(point, Sz); + point *Zptr = gZ; + double ZSum = 0.; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + if(p->r > 1.) continue; // throw out points with R>1 + double theta = p->theta; + // components of grad Zj: + + // 1. Theta_j; 2. dTheta_j/Dtheta + double Tj = 1., dTj = 0.; + if(m){ + double costm, sintm; + sincos(theta*(double)m_abs, &sintm, &costm); + if(m > 0){ + Tj = costm; + dTj = -m_abs * sintm; + }else{ + Tj = sintm; + dTj = m_abs * costm; + } + } + + // 3. R_j & dR_j/dr + double Rj = 0., dRj = 0.; + for(k = 0; k <= iup; k++){ + double rj = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//----------------------------------- ----- ------------------------------- + FK[k]*FK[isum-k]*FK[iup-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)! + ); +//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k); + int kidx = n - 2*k; + Rj += rj * Rpow[kidx][j]; // *R^{n-2k} + if(kidx > 0) + dRj += rj * kidx * Rpow[kidx - 1][j]; // *(n-2k)*R^{n-2k-1} + } + // normalisation for Zernike + double eps_m = (m) ? 1. : 2., sq = sqrt(2.*(double)(n+1) / M_PI / eps_m); + Rj *= sq, dRj *= sq; + // 4. sin/cos + double sint, cost; + sincos(theta, &sint, &cost); + + // projections of gradZj + double TdR = Tj * dRj, RdT = Rj * dTj / p->r; + Zptr->x = TdR * cost - RdT * sint; + Zptr->y = TdR * sint + RdT * cost; + // norm factor + ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y; + } + if(norm) *norm = ZSum; + // free unneeded memory + free_rpow(&Rpow, n); + return gZ; +} + +/** + * Build components of vector polynomial Sj + * @param p - index of polynomial + * @param Sz - number of points + * @param P (i) - coordinates of reference points + * @param norm (o) - norm factor or NULL + * @return Sj(n,m) on image points + */ +point *zerngradR(int p, int Sz, polar *P, double *norm){ + int n, m, i; + convert_Zidx(p, &n, &m); + if(n < 1) errx(1, "Order of gradient Z must be > 0!"); + int m_abs = iabs(m); + double Zsum, K = 1.; + point *Sj = gradZR(n, m, Sz, P, &Zsum); + if(n != m_abs && n > 2){ // avoid trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0 + K = sqrt(((double)n+1.)/(n-1.)); + Zsum = 0.; + point *Zj= gradZR(n-2, m, Sz, P, NULL); + point *Sptr = Sj, *Zptr = Zj; + for(i = 0; i < Sz; i++, Sptr++, Zptr++){ + Sptr->x -= K * Zptr->x; + Sptr->y -= K * Zptr->y; + Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y; + } + FREE(Zj); + } + if(norm) *norm = Zsum; + return Sj; +} + + +/** + * Decomposition of image with normals to wavefront by Zhao's polynomials + * by least squares method through LU-decomposition + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param grads(i) - wavefront gradients values in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of coefficients + */ +double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx){ + int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + if(Zsz) *Zsz = pmax; + // Now allocate matrix: its Nth row is equation for Nth data point, + // Mth column is Z_M + int Sz2 = Sz*2, x, y; + gsl_matrix *M = gsl_matrix_calloc(Sz2, pmax); + // fill matrix with coefficients + size_t T = M->tda, S = T * Sz; // T is matrix stride, S - index of second coefficient + for(x = 1; x < pmax; x++){ + double norm; + int n, m; + convert_Zidx(x, &n, &m); + point *dZcoeff = gradZR(n, m, Sz, P, &norm), *dZptr = dZcoeff; + //point *dZcoeff = zerngradR(x, Sz, P, &norm), *dZptr = dZcoeff; + double *ptr = &(M->data[x]); + // X-component is top part of resulting matrix, Y is bottom part + for(y = 0; y < Sz; y++, ptr+=T, dZptr++){ // fill xth polynomial + *ptr = dZptr->x; + ptr[S] = dZptr->y; + } + FREE(dZcoeff); + } + + gsl_vector *ans = gsl_vector_calloc(pmax); + gsl_vector *b = gsl_vector_calloc(Sz2); + double *ptr = b->data; + for(x = 0; x < Sz; x++, ptr++, grads++){ + // fill components of vector b like components of matrix M + *ptr = grads->x; + ptr[Sz] = grads->y; + } + + gsl_vector *tau = gsl_vector_calloc(pmax); + gsl_linalg_QR_decomp(M, tau); + + gsl_vector *residual = gsl_vector_calloc(Sz2); + gsl_linalg_QR_lssolve(M, tau, b, ans, residual); + ptr = &ans->data[1]; + int maxIdx = 0; + double *Zidxs = MALLOC(double, pmax); + for(x = 1; x < pmax; x++, ptr++){ + if(fabs(*ptr) < Z_prec) continue; + Zidxs[x] = *ptr; + maxIdx = x; + } + + gsl_matrix_free(M); + gsl_vector_free(ans); + gsl_vector_free(b); + gsl_vector_free(tau); + gsl_vector_free(residual); + + if(lastIdx) *lastIdx = maxIdx; + return Zidxs; +} + +/** + * Decomposition of image with normals to wavefront by Zhao's polynomials + * @param Nmax (i) - maximum power of Zernike polinomial for decomposition + * @param Sz, P(i) - size (Sz) of points array (P) + * @param grads(i) - wavefront gradients values in points P + * @param Zsz (o) - size of Z coefficients array + * @param lastIdx(o) - (if !NULL) last non-zero coefficient + * @return array of coefficients + */ +double *gradZdecomposeR(int Nmax, int Sz, polar *P, point *grads, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + double *Zidxs = MALLOC(double, pmax); + point *icopy = MALLOC(point, Sz); + memcpy(icopy, grads, Sz*sizeof(point)); // make image copy to leave it unchanged + *Zsz = pmax; + for(i = 1; i < pmax; i++){ // now we fill array + double norm; + point *dZcoeff = zerngradR(i, Sz, P, &norm); + int j; + point *iptr = icopy, *zptr = dZcoeff; + double K = 0.; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += zptr->x*iptr->x + zptr->y*iptr->y; // multiply matrixes to get coefficient + K /= norm; + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = dZcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values + iptr->y -= K * zptr->y; + } + FREE(dZcoeff); + } + if(lastIdx) *lastIdx = maxIdx; + FREE(icopy); + return Zidxs; +} + +/** + * Restoration of image with normals Zhao's polynomials coefficients + * all like Zcompose + * @return restored image + */ +point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + point *image = MALLOC(point, Sz); + for(i = 1; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + point *Zcoeff = zerngradR(i, Sz, P, NULL); + int j; + point *iptr = image, *zptr = Zcoeff; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + iptr->x += K * zptr->x; + iptr->y += K * zptr->y; + } + FREE(Zcoeff); + } + return image; +} diff --git a/Zernike/zernike_annular.c b/Zernike/zernike_annular.c new file mode 100644 index 0000000..10ff7ae --- /dev/null +++ b/Zernike/zernike_annular.c @@ -0,0 +1,386 @@ +/* + * zernike_annular.c - annular Zernike polynomials + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +/* + * These polynomials realized according to formulae in: +@ARTICLE{1981JOSA...71...75M, + author = {{Mahajan}, V.~N.}, + title = "{Zernike annular polynomials for imaging systems with annular pupils.}", + journal = {Journal of the Optical Society of America (1917-1983)}, + keywords = {Astronomical Optics}, + year = 1981, + volume = 71, + pages = {75-85}, + adsurl = {http://adsabs.harvard.edu/abs/1981JOSA...71...75M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@ARTICLE{2006ApOpt..45.3442H, + author = {{Hou}, X. and {Wu}, F. and {Yang}, L. and {Wu}, S. and {Chen}, Q. + }, + title = "{Full-aperture wavefront reconstruction from annular subaperture interferometric data by use of Zernike annular polynomials and a matrix method for testing large aspheric surfaces}", + journal = {\ao}, + keywords = {Mathematics, Optical data processing, Metrology, Optical testing}, + year = 2006, + month = may, + volume = 45, + pages = {3442-3455}, + doi = {10.1364/AO.45.003442}, + adsurl = {http://adsabs.harvard.edu/abs/2006ApOpt..45.3442H}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + */ + +#include "zernike.h" +#include "zern_private.h" + +// epsilon (==rmin) +double eps = -1., eps2 = -1.; // this global variable changes with each conv_r call! + +/** + * convert coordinates of points on ring into normalized coordinates in [0, 1] + * for further calculations of annular Zernike (expression for R_k^0) + * @param r0 (i) - original coordinates array + * @param Sz - size of r0 + * @return dinamically allocated array with new coordinates, which zero element is + * (0,0) + */ +polar *conv_r(polar *r0, int Sz){ + int x; + polar *Pn = MALLOC(polar, Sz); + double rmax = -1., rmin = 2.; + polar *p = r0, *p1 = Pn; + for(x = 0; x < Sz; x++, p++){ + if(p->r < rmin) rmin = p->r; + if(p->r > rmax) rmax = p->r; + } + if(rmin > rmax || fabs(rmin - rmax) <= DBL_EPSILON || rmax > 1. || rmin < 0.) + errx(1, "polar coordinates should be in [0,1]!"); + eps = rmin; + eps2 = eps*eps; + p = r0; + //rmax = 1. - eps2; // 1 - eps^2 -- denominator + rmax = rmax*rmax - eps2; + for(x = 0; x < Sz; x++, p++, p1++){ + p1->theta = p->theta; + p1->r = sqrt((p->r*p->r - eps2) / rmax); // r = (r^2 - eps^2) / (1 - eps^2) + } + return Pn; +} + +/** + * Allocate 2D array (H arrays of length W) + * @param len - array's width (number of elements in each 1D array + * @param num - array's height (number of 1D arrays) + * @return dinamically allocated 2D array filled with zeros + */ +_2D *_2ALLOC(int len, int num){ + _2D *r = MALLOC(_2D, 1); + r->data = MALLOC(double*, num); + int x; + for(x = 0; x < num; x++) + r->data[x] = MALLOC(double, len); + r->len = (size_t)len; + r->num = (size_t)num; + return r; +} + +void _2FREE(_2D **a){ + int x, W = (*a)->num; + for(x = 0; x < W; x++) free((*a)->data[x]); + free((*a)->data); + FREE(*a); +} + +/** + * Zernike radial function R_n^0 + * @param n - order of polynomial + * @param Sz - number of points + * @param P (i) - modified polar coordinates: sqrt((r^2-eps^2)/(1-eps^2)) + * @param R (o) - array with size Sz to which polynomial will be stored + */ +void Rn0(int n, int Sz, polar *P, double *R){ + if(Sz < 1 || !P) + errx(1, "Size of matrix must be > 0!"); + if(n < 0) + errx(1, "Order of Zernike polynomial must be > 0!"); + if(!FK) build_factorial(); + int j, k, iup = n/2; + double **Rpow = build_rpowR(n, Sz, P); + // now fill output matrix + double *Zptr = R; + polar *p = P; + for(j = 0; j < Sz; j++, p++, Zptr++){ + double Z = 0.; + if(p->r > 1.) continue; // throw out points with R>1 + // calculate R_n^m + for(k = 0; k <= iup; k++){ // Sum + double F = FK[iup-k]; + double z = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)! + /(//------------------------------------- --------------------- + FK[k]*F*F // k!((n/2-k)!)^2 + ); + Z += z * Rpow[n-2*k][j]; // *R^{n-2k} + } + *Zptr = Z; +//DBG("order %d, point (%g, %g): %g\n", n, p->r, p->theta, Z); + } + // free unneeded memory + free_rpow(&Rpow, n); +} + +/** + * Convert parameters n, m into index in Noll notation, p + * @param n, m - Zernike orders + * @return order in Noll notation + */ +inline int p_from_nm(int n, int m){ + return (n*n + 2*n + m)/2; +} + +/** + * Create array [p][pt] of annular Zernike polynomials + * @param pmax - max number of polynomial in Noll notation + * @param Sz - size of points array + * @param P (i) - array with points' coordinates + * @param Norm(o)- array [0..pmax] of sum^2 for each polynomial (dynamically allocated) or NULL + * @return + */ +#undef DBG +#define DBG(...) +_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm){ + int m, n, j, nmax, mmax, jmax, mW, jm; + convert_Zidx(pmax, &nmax, &mmax); + mmax = nmax; // for a case when p doesn't include a whole n powers + jmax = nmax / 2; + mW = mmax + 1; // width of row + jm = jmax + mW; // height of Q0 and h + polar *Pn = conv_r(P, Sz); + if(eps < 0. || eps2 < 0. || eps > 1. || eps2 > 1.) + errx(1, "Wrong epsilon value!!!"); + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + _2D *h = _2ALLOC(mW, jm); // constants -> array [(jmax+1)+((mmax+1)-1)][mmax+1] + _2D *Q0= _2ALLOC(mW, jm); // Q(0) -> array [(jmax+1)+((mmax+1)-1)][mmax+1] + _2D *Q = _2ALLOC(Sz, (jmax + 1)*mW); // functions -> array [(jmax+1)*(mmax+1)][Sz] +DBG("Allocated Q&h, size of Q: %zdx%zd, size of H: %zdx%zd\n", Q->num, Q->len, h->num, h->len); +DBG("pmax = %d, nmax = mmax = %d, jmax = %d, jm = %d\n", pmax, nmax, jmax, jm); +DBG("eps=%g, eps2=%g\n", eps, eps2); + // in Q index of [j][m][point] is [j*mmax+m][point] + double **hd = h->data, **Qd = Q->data, **Q0d = Q0->data; + // fill h_j^0 & Q_j^0 + for(j = 0; j < jm; j++){ // main matrix + triangle +DBG("j = %d, m = 0\n", j); + polar Zero = {0., 0.}; + // h_j^0 = (1-eps^2)/(2(2j+1)) + hd[j][0] = (1.-eps2)/(2.*j+1.)/2.; + // Q_j^0(rho, eps) = R_{2j}^0(r) +DBG("Q,h [%d][0], h=%g\n", j, hd[j][0]); + if(j <= jmax) Rn0(2*j, Sz, Pn, Qd[mW*j]); + Rn0(2*j, 1, &Zero, &Q0d[j][0]); // fill also Q(0) + } + FREE(Pn); + // fill h_j^m & Q_j&m + int JMAX = jm-1; // max number of j + for(m = 1; m <= mmax; m++, JMAX--){ // m in outern cycle because of Q_{j+1}! +DBG("\n\nm=%d\n",m); + int mminusone = m - 1; + int idx = mminusone;// index Q_j^{m-1} + for(j = 0; j < JMAX; j++, idx+=mW){ +DBG("\nj=%d\n",j); + // 2(2j+2m-1) * h_j^{m-1} + // H_j^m = ------------------------------ + // (j+m)(1-eps^2) * Q_j^{m-1}(0) + double H = 2*(2*(j+m)-1)/(j+m)/(1-eps2) * hd[j][mminusone]/Q0d[j][mminusone]; + // h_j^m = -H_j^m * Q_{j+1}^{m-1}(0) + hd[j][m] = -H * Q0d[j+1][mminusone]; +DBG("hd[%d][%d] = %g (H = %g)\n", j, m, hd[j][m], H); + // j Q_i^{m-1}(0) * Q_i^{m-1}(r) + // Q_j^m(r) = H_j^m * Sum ----------------------------- + // i=0 h_i^{m-1} + if(j <= jmax){ + int pt; + for(pt = 0; pt < Sz; pt++){ // cycle by points + int i, idx1 = mminusone; + double S = 0; + for(i = 0; i <=j; i++, idx1+=mW){ // cycle Sum + S += Q0d[i][mminusone] * Qd[idx1][pt] / hd[i][mminusone]; + } + Qd[idx+1][pt] = H * S; // Q_j^m(rho^2) +//DBG("Q[%d][%d](%d) = %g\n", j, m, pt, H * S); + } + } + // and Q(0): + double S = 0.; + int i; + for(i = 0; i <=j; i++){ + double qim = Q0d[i][mminusone]; + S += qim*qim / hd[i][mminusone]; + } + Q0d[j][m] = H * S; // Q_j^m(0) + DBG("Q[%d][%d](0) = %g\n", j, m, H * S); + } + } + _2FREE(&Q0); + + // allocate memory for Z_p + _2D *Zarr = _2ALLOC(Sz, pmax + 1); // pmax is max number _including_ + double **Zd = Zarr->data; + // now we should calculate R_n^m +DBG("R_{2j}^0\n"); + // R_{2j}^0(rho) = Q_j^0(rho^2) + size_t S = Sz*sizeof(double); + for(j = 0; j <= jmax; j++){ + int pidx = p_from_nm(2*j, 0); + if(pidx > pmax) break; + DBG("R[%d]\n", pidx); + memcpy(Zd[pidx], Qd[j*mW], S); + } +DBG("R_n^n\n"); + // R_n^n(rho) = rho^n / sqrt(Sum_{i=0}^n[eps^{2i}]) + double **Rpow = build_rpowR(nmax, Sz, P); // powers of rho + double epssum = 1., epspow = 1.; + for(n = 1; n <= nmax; n++){ + int pidx = p_from_nm(n, -n); // R_n^{-n} + if(pidx > pmax) break; + DBG("R[%d] (%d, %d)\n", pidx, n, -n); + epspow *= eps2; // eps^{2n} + epssum += epspow;// sum_0^n eps^{2n} + double Sq = sqrt(epssum); // sqrt(sum...) + double *rptr = Zd[pidx], *pptr = Rpow[n]; + int pt; + for(pt = 0; pt < Sz; pt++, rptr++, pptr++) + *rptr = *pptr / Sq; // rho^n / sqrt(...) + int pidx1 = p_from_nm(n, n); // R_n^{n} + if(pidx1 > pmax) break; + DBG("R[%d] (%d, %d)\n", pidx1, n, n); + memcpy(Zd[pidx1], Zd[pidx], S); + } +DBG("R_n^m\n"); + /* / 1 - eps^2 \ + R_{2j+|m|}^{|m|}(rho, eps) = sqrt| ------------------ | * rho^m * Q_j^m(rho^2) + \ 2(2j+|m|+1)*h_j^m / */ + for(j = 1; j <= jmax; j++){ + for(m = 1; m <= mmax; m++){ + int N = 2*j + m, pt, pidx = p_from_nm(N, -m); + if(pidx > pmax) continue; + DBG("R[%d]\n", pidx); + double *rptr = Zd[pidx], *pptr = Rpow[m], *qptr = Qd[j*mW+m]; + double hjm = hd[j][m]; + for(pt = 0; pt < Sz; pt++, rptr++, pptr++, qptr++) + *rptr = sqrt((1-eps2)/2/(N+1)/hjm) * (*pptr) * (*qptr); + int pidx1 = p_from_nm(N, m); + if(pidx1 > pmax) continue; + DBG("R[%d]\n", pidx1); + memcpy(Zd[pidx1], Zd[pidx], S); + } + } + free_rpow(&Rpow, nmax); + _2FREE(&h); + _2FREE(&Q); + // last step: fill Z_n^m and find sum^2 + int p; + double *sum2 = MALLOC(double, pmax+1); + for(p = 0; p <= pmax; p++){ + convert_Zidx(p, &n, &m); + DBG("p=%d (n=%d, m=%d)\n", p, n, m); + int pt; + double (*fn)(double val); + double empty(double val){return 1.;} + if(m == 0) fn = empty; + else{if(m > 0) + fn = cos; + else + fn = sin;} + double m_abs = fabs((double)m), *rptr = Zd[p], *sptr = &sum2[p]; + polar *Pptr = P; + for(pt = 0; pt < Sz; pt++, rptr++, Pptr++){ + *rptr *= fn(m_abs * Pptr->theta); + //DBG("\tpt %d, val = %g\n", pt, *rptr); + *sptr += (*rptr) * (*rptr); + } + DBG("SUM p^2 = %g\n", *sptr); + } + if(Norm) *Norm = sum2; + else FREE(sum2); + return Zarr; +} +#undef DBG +#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0) + +/** + * Compose wavefront by koefficients of annular Zernike polynomials + * @params like for ZcomposeR + * @return composed wavefront (dynamically allocated) + */ +double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){ + int i; + double *image = MALLOC(double, Sz); + _2D *Zannular = ann_Z(Zsz-1, Sz, P, NULL); + double **Zcoeff = Zannular->data; + for(i = 0; i < Zsz; i++){ // now we fill array + double K = Zidxs[i]; + if(fabs(K) < Z_prec) continue; + int j; + double *iptr = image, *zptr = Zcoeff[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++){ + *iptr += K * (*zptr); // add next Zernike polynomial + } + } + _2FREE(&Zannular); + return image; +} + + +/** + * Decompose wavefront by annular Zernike + * @params like ZdecomposeR + * @return array of Zernike coefficients + */ +double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){ + int i, pmax, maxIdx = 0; + double *Zidxs = NULL, *icopy = NULL; + pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation + Zidxs = MALLOC(double, pmax); + icopy = MALLOC(double, Sz); + memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged + *Zsz = pmax; + double *norm; + _2D *Zannular = ann_Z(pmax-1, Sz, P, &norm); + double **Zcoeff = Zannular->data; + for(i = 0; i < pmax; i++){ // now we fill array + int j; + double *iptr = icopy, *zptr = Zcoeff[i], K = 0., norm_i = norm[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++) + K += (*zptr) * (*iptr) / norm_i; // multiply matrixes to get coefficient + if(fabs(K) < Z_prec) + continue; // there's no need to substract values that are less than our precision + Zidxs[i] = K; + maxIdx = i; + iptr = icopy; zptr = Zcoeff[i]; + for(j = 0; j < Sz; j++, iptr++, zptr++) + *iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values + } + if(lastIdx) *lastIdx = maxIdx; + _2FREE(&Zannular); + FREE(norm); + FREE(icopy); + return Zidxs; +} diff --git a/bidirectional_list.c b/bidirectional_list.c new file mode 100644 index 0000000..227a473 --- /dev/null +++ b/bidirectional_list.c @@ -0,0 +1,139 @@ +/* + * bidirectional_list.c - bidi sorted list + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include + +#include "bidirectional_list.h" + +#ifndef Malloc +#define Malloc my_alloc +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} +#endif +#ifndef FREE +#define FREE(arg) do{free(arg); arg = NULL;}while(0) +#endif + + +List *l_search(List *root, int v, List **Left, List **Right){ + List *L = NULL, *R = NULL; + int dir = -1; + if(!root) goto not_Fnd; + if(root->data == v) + return root; + if(root->data < v) dir = 1; + do{ + L = root->left_p; // left leaf + R = root->right_p; // right leaf + int D = root->data;// data in leaf + if(D > v && dir == -1){ // search to left is true + R = root; + root = L; + } + else if(D < v && dir == 1){ // search to right is true + L = root; + root = R; + } + else if(D == v){ // we found exact value + if(Left) *Left = L; + if(Right) *Right = R; + return root; + } + else break; + }while(root); // if root == NULL, we catch an edge + // exact value not found + if(root){ // we aren't on an edge + if(dir == -1) L = root; + else R = root; + } +not_Fnd: + if(Left) *Left = L; + if(Right) *Right = R; + return NULL; // value not found +} + + +List *l_insert(List *root, int v){ + List *node, *L, *R; + node = l_search(root, v, &L, &R); + if(node) return node; // we found exact value V + // there's no exact value: insert new + if((node = (List*) Malloc(1, sizeof(List))) == 0) return NULL; // allocation error + node->data = v; // insert data + // add neighbours: + node->left_p = L; + node->right_p = R; + // fix neighbours: + if(L) L->right_p = node; + if(R) R->left_p = node; + return node; +} + +void l_remove(List **root, int v){ + List *node, *left, *right; + if(!root) return; + if(!*root) return; + node = l_search(*root, v, &left, &right); + if(!node) return; // there's no such element + if(node == *root) *root = node->right_p; + if(!*root) *root = node->left_p; + FREE(node); + if(left) left->right_p = right; + if(right) right->left_p = left; +} + +#ifdef STANDALONE +int main(void) { + List *tp, *root_p = NULL; + int i, ins[] = {4,2,6,1,3,4,7}; + for(i = 0; i < 7; i++) + if(!(root_p = l_insert(root_p, ins[i]))) + err(1, "Malloc error"); // can't insert + for(i = 0; i < 9; i++){ + tp = l_search(root_p, i, NULL, NULL); + printf("%d ", i); + if(!tp) printf("not "); + printf("found\n"); + } + printf("now delete items 1, 4 and 8\n"); + l_remove(&root_p, 1); + l_remove(&root_p, 4); + l_remove(&root_p, 8); + for(i = 0; i < 9; i++){ + tp = l_search(root_p, i, NULL, NULL); + printf("%d ", i); + if(!tp) printf("not "); + printf("found\n"); + } + return 0; +} +#endif diff --git a/bidirectional_list.h b/bidirectional_list.h new file mode 100644 index 0000000..28b8379 --- /dev/null +++ b/bidirectional_list.h @@ -0,0 +1,34 @@ +/* + * bidirectional_list.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __BIDIRECTIONAL_LIST_H__ +#define __BIDIRECTIONAL_LIST_H__ + +typedef struct list_node{ + int data; + struct list_node *left_p, *right_p; +} List; + +List *l_search(List *root, int v, List **Left, List **Right); +List *l_insert(List *root, int v); +void l_remove(List **root, int v); + +#endif // __BIDIRECTIONAL_LIST_H__ diff --git a/cmakelists_/CMakeLists_regular_01.txt b/cmakelists_/CMakeLists_regular_01.txt new file mode 100644 index 0000000..fbe7687 --- /dev/null +++ b/cmakelists_/CMakeLists_regular_01.txt @@ -0,0 +1,111 @@ +cmake_minimum_required(VERSION 2.8) +set(PROJ tvguide) +set(MINOR_VERSION "0") +set(MID_VERSION "0") +set(MAJOR_VERSION "0") +set(VERSION "${MAJOR_VERSION}.${MID_VERSION}.${MINOR_VERSION}") + +message("VER: ${VERSION}") + +# default flags +set(CFLAGS -O2 -Wextra -Wall -Werror -W -std=gnu99) + +set(CMAKE_COLOR_MAKEFILE ON) + +# here is one of two variants: all .c in directory or .c files in list +aux_source_directory(${CMAKE_CURRENT_SOURCE_DIR} SOURCES) +#set(SOURCES list_of_c_files) + +# we can change file list +#if(NOT DEFINED something) +# set(SOURCES ${SOURCES} one_more_list) +# add_definitions(-DSOME_DEFS) +#endif() + +# cmake -DDEBUG=1 -> debugging +if(DEFINED DEBUG) + add_definitions(-DEBUG) +endif() + +# directory should contain dir locale/ru for gettext translations +set(LCPATH ${CMAKE_SOURCE_DIR}/locale/ru) + +if(NOT DEFINED LOCALEDIR) + if(DEFINED DEBUG) + set(LOCALEDIR ${CMAKE_CURRENT_SOURCE_DIR}/locale) + else() + set(LOCALEDIR ${CMAKE_INSTALL_PREFIX}/share/locale) + endif() +endif() + +# gettext modules +set(MODULES libavformat libavcodec libswscale libavutil libavdevice) +# additional modules on condition +#if(DEFINED SOMETHING) +# set(MODULES ${MODULES} more_modules>=version) +# add_definitions(-DSOMEDEFS) +#endif() + +project(${PROJ}) +# change wrong behaviour with install prefix +if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT AND CMAKE_INSTALL_PREFIX MATCHES "/usr/local") + message("Change default install path to /usr") + set(CMAKE_INSTALL_PREFIX "/usr") +endif() +message("Install dir prefix: ${CMAKE_INSTALL_PREFIX}") + +# gettext files +set(PO_FILE ${LCPATH}/messages.po) +set(MO_FILE ${LCPATH}/LC_MESSAGES/${PROJ}.mo) +set(RU_FILE ${LCPATH}/ru.po) + +# pkgconfig +find_package(PkgConfig REQUIRED) +pkg_check_modules(${PROJ} REQUIRED ${MODULES}) + +# exe file +add_executable(${PROJ} ${SOURCES} ${PO_FILE} ${MO_FILE}) +target_link_libraries(${PROJ} ${${PROJ}_LIBRARIES}) +include_directories(${${PROJ}_INCLUDE_DIRS}) +link_directories(${${PROJ}_LIBRARY_DIRS}) +add_definitions(${CFLAGS} -DLOCALEDIR=\"${LOCALEDIR}\" + -DPACKAGE_VERSION=\"${VERSION}\" -DGETTEXT_PACKAGE=\"${PROJ}\" + -DMINOR_VERSION=\"${MINOR_VERSION}\" -DMID_VERSION=\"${MID_VERSION}\" + -DMAJOR_VERSION=\"${MAJOR_VESION}\") + +# Installation of the program +INSTALL(FILES ${MO_FILE} DESTINATION "share/locale/ru/LC_MESSAGES") + #PERMISSIONS OWNER_WRITE OWNER_READ GROUP_READ WORLD_READ) +INSTALL(TARGETS ${PROJ} DESTINATION "bin") + #PERMISSIONS OWNER_WRITE OWNER_READ OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE) +# Script to be executed at installation time (kind of post-intallation script) to +# change the right accesses on the installed files +#INSTALL(SCRIPT inst.cmake) + +find_package(Gettext REQUIRED) +find_program(GETTEXT_XGETTEXT_EXECUTABLE xgettext) +if(NOT GETTEXT_XGETTEXT_EXECUTABLE OR NOT GETTEXT_MSGFMT_EXECUTABLE) + message(FATAL_ERROR "xgettext not found") +endif() +file(MAKE_DIRECTORY ${LCPATH}) +file(MAKE_DIRECTORY ${LCPATH}/LC_MESSAGES) + +add_custom_command( + OUTPUT ${PO_FILE} + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + COMMAND ${GETTEXT_XGETTEXT_EXECUTABLE} --from-code=utf-8 ${SOURCES} -c -k_ -kN_ -o ${PO_FILE} + COMMAND sed -i 's/charset=.*\\\\n/charset=koi8-r\\\\n/' ${PO_FILE} + COMMAND enconv ${PO_FILE} + DEPENDS ${SOURCES} +) +# we need this to prewent ru.po from deleting by make clean +add_custom_target( + RU_FILE + COMMAND [ -f ${RU_FILE} ] && ${GETTEXT_MSGMERGE_EXECUTABLE} -Uis ${RU_FILE} ${PO_FILE} || cp ${PO_FILE} ${RU_FILE} + DEPENDS ${PO_FILE} +) +add_custom_command( + OUTPUT ${MO_FILE} + COMMAND make RU_FILE && ${GETTEXT_MSGFMT_EXECUTABLE} ${RU_FILE} -o ${MO_FILE} + DEPENDS ${PO_FILE} +) diff --git a/erosion-dilation/Makefile b/erosion-dilation/Makefile new file mode 100644 index 0000000..15f4766 --- /dev/null +++ b/erosion-dilation/Makefile @@ -0,0 +1,19 @@ +PROGRAM = binmorph_test +CXX.SRCS = binmorph.c main.c +CC = gcc +LDFLAGS = -lm -fopenmp +CXX = gcc +DEFINES = -DEBUG -DOMP +CFLAGS = -fopenmp -Wall -Werror -O3 -std=gnu99 $(DEFINES) +OBJS = $(CXX.SRCS:.c=.o) +all : $(PROGRAM) clean +binmorph.c : cclabling.h cclabling_1.h dilation.h fs_filter.h fc_filter.h +$(PROGRAM) : $(OBJS) + $(CC) $(LDFLAGS) $(OBJS) -o $(PROGRAM) +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) + +### +# name1.o : header1.h header2.h ... diff --git a/erosion-dilation/README b/erosion-dilation/README new file mode 100644 index 0000000..90cde95 --- /dev/null +++ b/erosion-dilation/README @@ -0,0 +1,523 @@ +************************************ + / +************************************ + + + + +00000 00000 +01110 00000 +01110 => 00100 +01110 00000 +00000 00000 + + , 4 . +.. & . 4- : + +ABC +DEF +GHK + + E F, B H E. + E D, B E E. + + , : +1. E [E] ( ) +2. &: E = E & B & H +3. E D E F. + +, : (A&0X80)>>8 ( ) + + + + +00000 00110 +00110 01111 +01010 => 11111 +01100 11110 +00000 01100 + + : & | : +1. E ( ) +2. |: E = E | B | H +3. |. + + + + +************************************************** + +************************************************** + + : E. +ABC +DEF +GHK + + 4- + - , 4-, : + +E4 = E & ( B|H|(E<<1)|(E>>1)|(D<<7)|(F>>7) ) + +.. E , +4- . + +, ( /) E4 + ( - / B,H; + - / D, F). + + boost, - + , . + +, A,C,G K, +( ). + + + 8- : + + {x} = (x|(x<<1)|(x>>1)), +E8 = E & ( {B|H} | ((E<<1)|(E>>1)) | ((C|F|K)>>7) | ((A|D|G)<<7) ) + + , + . , , + {x} ( (x<<7) (x>>7) ???). +/* () , x<<7) (x>>7) + * E4 */ + + + , " " + , . , + 8- 8-. + +, - : + , . + , + ( , +, ). + + - () + (!). + - , , + , + . + + , + , ( "" +). + + +ABC +DEF +GHK + + , , + -> ( + , +). , +, . + + +=== === + +/* + , .. . +*/ + +, ( ) + . + E, 1. + , . , A,B,C D (.. - + ). F, G, H K. + E. - 1, + .. , , , + , .. , +- . + + , , , - +, . : + . . + +, , , + ! + + 1 , "" + ( 0). . + + : , + , , + ( - + ); , , + , "" . , : + , . + + +=== "". === + (, +, , , , + ). + +== == + , , , + : + +#include + +int incr(int a){ + fprintf(stderr, "%d ", a++); + incr(a); + return a; +} + +main(){ + incr(1); + return 0; +} + + + . : + +... 261504 (core dumped) +... 261350 (core dumped) +... 261333 (core dumped) + +, - . , , - . + 256 (262144): , ( + stdio). + + : 1 : + +#include +#include +#include + +int incr(int a){ + char *x = malloc(1024); + memset(x, 1, 1024); + fprintf(stderr, "%d ", a++); + incr(a); + return a; +} + +main(){ + incr(1); + return 0; +} + + " ": + +... 174258 (core dumped) +... 174372 (core dumped) +... 174325 (core dumped) + + 85? 85? `char x = 2;` - + . . + , - : + , 170 ( ), + . +## + + : ( - + ) ( + - , "" ). + + "" , + , . +, , + . + + , , "", + . + + +### " " ### + - , + "": + "". "". + "", ( - +, ), "" , . + , , + . , , + , "" 1 (.. + 1, - 1 ..). + 1, "" + 1 + .. , "" +. , , "" - +, , "" + . + .: "" +" ", ( + + 1 , .. ), + : 1 , "" "". + , , "" + . + + - . , + - . + +=== === +, . . + ( ?) : + (1 - , 2 - + , 3 - , 4 - +). , , + - ( , 8- +), - , - + , . +, : + . + + , == N. + + : + , . , + false . == true, , + - . + - "", 6 ( 8- +) ( - .. +, - - + , , , , + ). , + , : - ( - + , DL DR , +, ). + + , , + , ( , 1/0; + , . 256, ): +(X - , U - , V - ) + + A B C D E F G H I J K L M N O P Q R +XXX 000 UUU 000 101 101 11X 000 010 000 000 000 000 000 000 000 000 000 +11X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X +XXX 000 000 VVV 000 VVV VVV 010 000 000 000 000 000 000 000 000 000 000 + +A - +B - +C - , , + +D - , B + +E - + , + +F - +G - F: + +H - , .. I +I - H + +, : +- , , + : , , + - +; ( ) +- !=0, , + +- !=0 , , + + + : +X Y, Y Z, +X Y. +, + , assoc[assoc[X]] == assoc[X], .. + . + + : , +- . ( , +.. 5 ), , , . + . + + , , : +, "X", + "Y" (.. assoc[X] != X) "Z". +Z Y, , Z, <= Ncur , + Ncur. + , + , +" ": 3. + +, , +: . , , , + (1, 3 - , , +4 - ), , , ! + + +== == +, : + , ( - "" + "" + ( , , - +). + +, , + : +1. , 1, "" + "" ( ) + ; +2. ("" + 4 ); +3. , .1. + + "" "" .1, + " ", "" . .., "" + , . +, : + . , "" + , : + . + +, : . + , + + , 0. + , . +, . , ! + + - " " + -, - : + - - , - (, +, ). + +, - " " + . , +( - 4-, , .. - + , 4- +), CPU openMP + , : , +, . + + +// - : + FIFO , LIFO ... + + +== == + , - "". "" +inline- . + "" . + , - . + + -, + - - - "" : + , + , . + +, : + 4- 4-. + + "" , + . . - + . + -, , ( ) +"" ( FCfilter), . + +, 8- , + FCfilter, , + 4- , .. , 4- + . , + 8- ( ). + + : + ( ) . + . + . , + , . + + , 4- , . + . +8- , + . , , , + , . + "" +, , +, : + . , , + . + +! !!! + + : ( 20002000): , +cc4 2 , cc8!!! ? cc8 + ... ! +: 13.3 cc4 7.8 cc8. , + 4! cc4 `FC_filter`, + (, , +, - , ). + + , , , + . cclabel. + , 2 4 , + . - ! + , . + - 200 ! + , 200000 +, , + ! + + , , . + . , , , + - , + , . , + - cc4 cc8! + + 50/50 10/90, + cc4: ! cc8 + 1. , : + ( ). + , : 50/50 + 20! 10/90 , ( + :( ). + + , - + ( ) + -> -. , +" " ( ). + + , , - . + - 11-. , 4:30 . + +== == + , . , + : () + , - . + , , : +- 1 2 : ; +- ( - ) , - + . +, , : + ( , ), + - , + . , + , 1/6 , " " + 1/4 . + + : , + , , , + : , + ( , , .. + ) , + , , - + . , , + . , + , - , + . + + () : + . , . + + - 8- , + 4- , . + + , - "", . , + - : ! + - . + , . + + : wifi , - , +! , : tracepath + ! , . + +- - . , + : + ! + +! ! ... diff --git a/erosion-dilation/binmorph.c b/erosion-dilation/binmorph.c new file mode 100644 index 0000000..eb6b8f3 --- /dev/null +++ b/erosion-dilation/binmorph.c @@ -0,0 +1,428 @@ +/* + * binmorph.c - functions for morphological operations on binary image + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + + +/* + TEST (10000 dilations / erosions / substitutions) for i5-3210M CPU @ 2.50GHz x2cores x2hyper + 1. no openmp: + * 28.8 / 29.2 / 14.8 + 2. 2 threads: + * 17.0 / 17.9 / 8.8 + 3. 4 threads: + * 17.0 / 17.5 / 8.8 + 4. 8 threads: + * 17.0 / 17.7 / 8.9 + * + option -O3 gives for nthreads = 4: + * 6.9 / 6.9 / 2.9 +*/ + +#include +#include +#include +#include +#include +#include "binmorph.h" +//#include "fifo.h" + +#ifdef OMP + #ifndef OMP_NUM_THREADS + #define OMP_NUM_THREADS 4 + #endif + #define Stringify(x) #x + #define OMP_FOR(x) _Pragma(Stringify(omp parallel for x)) +#else + #define OMP_FOR(x) +#endif + +// global arrays for erosion/dilation masks +unsigned char *ER = NULL, *DIL = NULL; +bool __Init_done = false; // == true if arrays are inited + +/* + * =================== AUXILIARY FUNCTIONS ===================> + */ + +/** + * function for different purposes that need to know time intervals + * @return double value: time in seconds + */ +double dtime(){ + double t; + struct timeval tv; + gettimeofday(&tv, NULL); + t = tv.tv_sec + ((double)tv.tv_usec)/1e6; + return t; +} + +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} + +/** + * This function inits masks arrays for erosion and dilation + * You may call it yourself or it will be called when one of + * `erosion` or `dilation` functions will be ran first time + */ +void morph_init(){ + if(__Init_done) return; + int i;//,j; + //bool a[8], b[8], c[9]; + ER = Malloc(256, 1); + DIL = Malloc(256, 1); + for(i = 0; i < 256; i++){ + /*ER[i] = DIL[i] = 0; + for(j = 0; j < 8; j++){ + a[j] = (i >> j) & 1; + b[j] = a[j]; + c[j] = a[j]; + } + for(j = 1; j < 7; j++) + if(!a[j]) + b[j+1] = b[j-1]= false; + else + c[j] = c[j-1] = c[j+1] = true; + b[1] &= a[0]; b[6] &= a[7]; + c[1] |= a[0]; c[6] |= a[7]; + for(j = 0; j < 8; j++){ + ER[i] |= b[j] << j; + DIL[i] |= c[j] << j; + } + */ + ER[i] = i & ((i << 1) | 1) & ((i >> 1) | (0x80)); // don't forget that << and >> set borders to zero + DIL[i] = i | (i << 1) | (i >> 1); + } + __Init_done = true; +} + +/** + * Print out small "packed" images as matrix of "0" and "1" + * @param i (i) - image + * @param W, H - size of image + */ +void printC(unsigned char *i, int W, int H){ + int x,y,b; + for(y = 0; y < H; y++){ + for(x = 0; x < W; x++, i++){ + unsigned char c = *i; + for(b = 0; b < 8; b++) + printf("%c", c << b & 0x80 ? '1' : '0'); + } + printf("\n"); + } +} + +/** + * Print out small "unpacked" images as almost equiaxed matrix of ".." and "##" + * @param img (i) - image + * @param W, H - size of image + */ +void printB(bool *img, int W, int H){ + int i, j; + for(j = 0; j < W; j++) + if(j % 10 == 1){ printf("%02d", j-1); printf(" ");} + printf("\n"); + for(i = 0; i < H; i++){ + printf("%2d ",i); + for(j = 0; j < W; j++) + printf("%s", (*img++) ? "##" : ".."); + printf("\n"); + } +} + +/* + * <=================== AUXILIARY FUNCTIONS =================== + */ + +/* + * =================== CONVERT IMAGE TYPES ===================> + */ + +/** + * Convert boolean image into pseudo-packed (1 char == 8 pixels) + * @param im (i) - image to convert + * @param W, H - size of image im (must be larger than 1) + * @param W_0 (o) - (stride) new width of image + * @return allocated memory area with "packed" image + */ +unsigned char *bool2char(bool *im, int W, int H, int *W_0){ + if(W < 2 || H < 2) errx(1, "bool2char: image size too small"); + int y, W0 = (W + 7) / 8; + unsigned char *ret = Malloc(W0 * H, 1); + OMP_FOR() + for(y = 0; y < H; y++){ + int x, i, X; + bool *ptr = &im[y*W]; + unsigned char *rptr = &ret[y*W0]; + for(x = 0, X = 0; x < W0; x++, rptr++){ + for(i = 7; i > -1 && X < W; i--, X++, ptr++){ + *rptr |= *ptr << i; + } + } + } + if(W_0) *W_0 = W0; + return ret; +} + +/** + * Convert "packed" image into boolean + * @param image (i) - input image + * @param W, H, W_0 - size of image and width of "packed" image + * @return allocated memory area with "unpacked" image + */ +bool *char2bool(unsigned char *image, int W, int H, int W_0){ + int y; + bool *ret = Malloc(W * H, sizeof(bool)); + OMP_FOR() + for(y = 0; y < H; y++){ + int x, X, i; + bool *optr = &ret[y*W]; + unsigned char *iptr = &image[y*W_0]; + for(x = 0, X = 0; x < W_0; x++, iptr++) + for(i = 7; i > -1 && X < W; i--, X++, optr++){ + *optr = (*iptr >> i) & 1; + } + } + return ret; +} + +/** + * Convert "packed" image into size_t array for conncomp procedure + * @param image (i) - input image + * @param W, H, W_0 - size of image and width of "packed" image + * @return allocated memory area with copy of an image + */ +size_t *char2st(unsigned char *image, int W, int H, int W_0){ + int y; + size_t *ret = Malloc(W * H, sizeof(size_t)); + OMP_FOR() + for(y = 0; y < H; y++){ + int x, X, i; + size_t *optr = &ret[y*W]; + unsigned char *iptr = &image[y*W_0]; + for(x = 0, X = 0; x < W_0; x++, iptr++) + for(i = 7; i > -1 && X < W; i--, X++, optr++){ + *optr = (*iptr >> i) & 1; + } + } + return ret; +} + +/* + * <=================== CONVERT IMAGE TYPES =================== + */ + +/* + * =================== MORPHOLOGICAL OPERATIONS ===================> + */ + +/** + * Remove all non-4-connected pixels + * @param image (i) - input image + * @param W, H - size of image + * @return allocated memory area with converted input image + */ +unsigned char *FC_filter(unsigned char *image, int W, int H){ + if(W < 2 || H < 2) errx(1, "4-connect: image size too small"); + unsigned char *ret = Malloc(W*H, 1); + int y = 0, w = W-1, h = H-1; + // top of image, y = 0 + #define IM_UP + #include "fc_filter.h" + #undef IM_UP + // mid of image, y = 1..h-1 + #include "fc_filter.h" + // image bottom, y = h + y = h; + #define IM_DOWN + #include "fc_filter.h" + #undef IM_DOWN + return ret; +} + +/** + * Make morphological operation of dilation + * @param image (i) - input image + * @param W, H - size of image + * @return allocated memory area with dilation of input image + */ +unsigned char *dilation(unsigned char *image, int W, int H){ + if(W < 2 || H < 2) errx(1, "Dilation: image size too small"); + if(!__Init_done) morph_init(); + unsigned char *ret = Malloc(W*H, 1); + int y = 0, w = W-1, h = H-1; + // top of image, y = 0 + #define IM_UP + #include "dilation.h" + #undef IM_UP + // mid of image, y = 1..h-1 + #include "dilation.h" + // image bottom, y = h + y = h; + #define IM_DOWN + #include "dilation.h" + #undef IM_DOWN + return ret; +} + +/** + * Make morphological operation of erosion + * @param image (i) - input image + * @param W, H - size of image + * @return allocated memory area with erosion of input image + */ +unsigned char *erosion(unsigned char *image, int W, int H){ + if(W < 2 || H < 2) errx(1, "Erosion: image size too small"); + if(!__Init_done) morph_init(); + unsigned char *ret = Malloc(W*H, 1); + int y, w = W-1, h = H-1; + OMP_FOR() + for(y = 1; y < h; y++){ // reset first & last rows of image + unsigned char *iptr = &image[W*y]; + unsigned char *optr = &ret[W*y]; + unsigned char p = ER[*iptr] & 0x7f & iptr[-W] & iptr[W]; + int x; + if(!(iptr[1] & 0x80)) p &= 0xfe; + *optr++ = p; + iptr++; + for(x = 1; x < w; x++, iptr++, optr++){ + p = ER[*iptr] & iptr[-W] & iptr[W]; + if(!(iptr[-1] & 1)) p &= 0x7f; + if(!(iptr[1] & 0x80)) p &= 0xfe; + *optr = p; + } + p = ER[*iptr] & 0xfe & iptr[-W] & iptr[W]; + if(!(iptr[-1] & 1)) p &= 0x7f; + *optr++ = p; + iptr++; + } + return ret; +} + +/* + * <=================== MORPHOLOGICAL OPERATIONS =================== + */ + +/* + * =================== LOGICAL OPERATIONS ===================> + */ + +/** + * Logical AND of two images + * @param im1, im2 (i) - two images + * @param W, H - their size (of course, equal for both images) + * @return allocated memory area with image = (im1 AND im2) + */ +unsigned char *imand(unsigned char *im1, unsigned char *im2, int W, int H){ + unsigned char *ret = Malloc(W*H, 1); + int y; + OMP_FOR() + for(y = 0; y < H; y++){ + int x, S = y*W; + unsigned char *rptr = &ret[S], *p1 = &im1[S], *p2 = &im2[S]; + for(x = 0; x < W; x++) + *rptr++ = *p1++ & *p2++; + } + return ret; +} + +/** + * Substitute image 2 from image 1: reset to zero all bits of image 1 which set to 1 on image 2 + * @param im1, im2 (i) - two images + * @param W, H - their size (of course, equal for both images) + * @return allocated memory area with image = (im1 AND (!im2)) + */ +unsigned char *substim(unsigned char *im1, unsigned char *im2, int W, int H){ + unsigned char *ret = Malloc(W*H, 1); + int y; + OMP_FOR() + for(y = 0; y < H; y++){ + int x, S = y*W; + unsigned char *rptr = &ret[S], *p1 = &im1[S], *p2 = &im2[S]; + for(x = 0; x < W; x++) + *rptr++ = *p1++ & (~*p2++); + } + return ret; +} + +/* + * <=================== LOGICAL OPERATIONS =================== + */ + +/* + * =================== CONNECTED COMPONENTS LABELING ===================> + */ + + + +/* + * <=================== CONNECTED COMPONENTS LABELING =================== + */ +/* +#undef DBG +#undef EBUG +#define DBG(...) +*/ + +/** + * label 4-connected components on image + * (slow algorythm, but easy to parallel) + * + * @param I (i) - image ("packed") + * @param W,H,W_0 - size of the image (W - width in pixels, W_0 - width in octets) + * @param Nobj (o) - number of objects found + * @return an array of labeled components + */ +CCbox *cclabel4(unsigned char *Img, int W, int H, int W_0, size_t *Nobj){ + unsigned char *I = FC_filter(Img, W_0, H); + #include "cclabling.h" + FREE(I); + return ret; +} + +// label 8-connected components, look cclabel4 +CCbox *cclabel8(unsigned char *I, int W, int H, int W_0, size_t *Nobj){ + //unsigned char *I = EC_filter(Img, W_0, H); + #define LABEL_8 + #include "cclabling.h" + #undef LABEL_8 + // FREE(I); + return ret; +} + +CCbox *cclabel8_1(unsigned char *I, int W, int H, int W_0, size_t *Nobj){ + #define LABEL_8 + #include "cclabling_1.h" + #undef LABEL_8 + return ret; +} + +/* + * <=================== template ===================> + */ diff --git a/erosion-dilation/binmorph.h b/erosion-dilation/binmorph.h new file mode 100644 index 0000000..7761965 --- /dev/null +++ b/erosion-dilation/binmorph.h @@ -0,0 +1,86 @@ +/* + * binmorph.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __EROSION_DILATION_H__ +#define __EROSION_DILATION_H__ + +#include +#include + +#ifdef EBUG +#ifndef DBG + #define DBG(...) do{fprintf(stderr, __VA_ARGS__);}while(0) +#endif +#endif + +#define _U_ __attribute__((__unused__)) + +#ifndef Malloc +#define Malloc my_alloc +#endif + +#ifndef FREE +#define FREE(arg) do{free(arg); arg = NULL;}while(0) +#endif + +#ifndef _ +#define _(X) X +#endif + +// auxiliary functions +void *my_alloc(size_t N, size_t S); +double dtime(); +void printC(unsigned char *i, int W, int H); +void printB(bool *i, int W, int H); +void morph_init(); // there's no need to run it by hands, but your will is the law + +// convert image types +unsigned char *bool2char(bool *im, int W, int H, int *stride); +bool *char2bool(unsigned char *image, int W, int H, int W_0); +size_t *char2st(unsigned char *image, int W, int H, int W_0); + +// morphological operations +unsigned char *dilation(unsigned char *image, int W, int H); +unsigned char *erosion(unsigned char *image, int W, int H); +unsigned char *FC_filter(unsigned char *image, int W, int H); + +// logical operations +unsigned char *imand(unsigned char *im1, unsigned char *im2, int W, int H); +unsigned char *substim(unsigned char *im1, unsigned char *im2, int W, int H); + +// conncomp +// this is a box structure containing one object; data is aligned by original image bytes! +typedef struct { + unsigned char *data; // pattern data in "packed" format + int x, // x coordinate of LU-pixel of box in "unpacked" image (by pixels) + y, // y -//- + x_0; // x coordinate in "packed" image (morph operations should work with it) + size_t N;// number of component, starting from 1 +} CCbox; + +CCbox *cclabel4(unsigned char *I, int W, int H, int W_0, size_t *Nobj); +CCbox *cclabel8(unsigned char *I, int W, int H, int W_0, size_t *Nobj); + +CCbox *cclabel8_1(unsigned char *I, int W, int H, int W_0, size_t *Nobj); +#endif // __EROSION_DILATION_H__ + + diff --git a/erosion-dilation/cclabling.h b/erosion-dilation/cclabling.h new file mode 100644 index 0000000..744831b --- /dev/null +++ b/erosion-dilation/cclabling.h @@ -0,0 +1,204 @@ +/* + * cclabling.h - inner part of functions cclabel4 and cclabel8 + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +double t0 = dtime(); + CCbox *ret = NULL; + size_t N _U_ = 0, Ncur = 0; + size_t *labels = char2st(I, W, H, W_0); + int y; +#ifdef EBUG +for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + if(*ptr) printf("%02zx", *ptr); + else printf(" "); + } + printf("\n"); +} +FREE(labels); return ret; +#endif +printf("time for char2st: %gs\n", dtime()-t0); +t0 = dtime(); + // Step 1: mark neighbours by strings + size_t *ptr = labels; + for(y = 0; y < H; y++){ + bool found = false; + for(int x = 0; x < W; x++, ptr++){ + if(!*ptr){ found = false; continue;} + if(!found){ + found = true; + ++Ncur; + } + *ptr = Ncur; + } + } +printf("\n\ntime for step1: %gs (Ncur=%zd)\n", dtime()-t0, Ncur); +t0 = dtime(); + DBG("Initial mark\n"); + #ifdef EBUG +for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + if(*ptr) printf("%02zx", *ptr); + else printf(" "); + } + printf("\n"); +} + #endif + + // Step 2: fill associative array with remarking + DBG("remarking\n"); + N = Ncur+1; // size of markers array (starting from 1) + size_t i, *assoc = Malloc(N, sizeof(size_t)); + for(i = 0; i < N; i++) assoc[i] = i; // primary initialisation + // now we should scan image again to rebuild enumeration + // THIS PROCESS AVOID PARALLELISATION??? + Ncur = 0; + size_t h = H-1 + #ifdef LABEL_8 + ,w = W-1; + #endif + ; + inline void remark(size_t old, size_t *newv){ + size_t new = *newv; + if(assoc[old] == new){ + DBG("(%zd)\n", new); + return; + } + DBG("[%zx], %zx -> %zx ", Ncur, old, new); + if(assoc[old] == old){ // remark non-remarked value + assoc[old] = new; + DBG("\n"); + return; + } + // value was remarked -> we have to remark "new" to assoc[old] + // and decrement all marks, that are greater than "new" + size_t ao = assoc[old]; + DBG("(remarked to %zx) ", assoc[old]); + // now we must check assoc[old]: if it is < new -> change *newv, else remark assoc[old] + if(ao < new) + *newv = ao; + else{ + size_t x = ao; ao = new; new = x; // swap values + } + int xx = old + W / 2; + if(xx > N) xx = N; + OMP_FOR() + for(int i = 1; i <= xx; i++){ + size_t m = assoc[i]; + if(m < new) continue; + if(m == new){ + assoc[i] = ao; + DBG("remark %x (%zd) to %zx ", i, m, ao); + }else if(m <= Ncur){ + assoc[i]--; + DBG("decr %x: %zx, ", i, assoc[i]); + } + } + DBG("\n"); + Ncur--; + } + for(y = 0; y < H; y++){ // FIXME!!!! throw out that fucking "if" checking coordinates!!! + size_t *ptr = &labels[y*W]; + bool found = false; + for(int x = 0; x < W; x++, ptr++){ + size_t curval = *ptr; + if(!curval){ found = false; continue;} + if(found) continue; // we go through remarked pixels + found = true; + DBG("curval: %zx ", curval); + // now check neighbours in upper and lower string: + size_t *U = (y) ? &ptr[-W] : NULL; + size_t *D = (y < h) ? &ptr[W] : NULL; + size_t upmark = 0; // mark from line above + if(U){ + #ifdef LABEL_8 + if(x && U[-1]){ // there is something in upper left corner -> make remark by its value + upmark = assoc[U[-1]]; + }else // check point above only if there's nothing in left up + #endif + if(U[0]) upmark = assoc[U[0]]; + #ifdef LABEL_8 + if(x < w) if(U[1]){ // there's something in upper right + if(upmark){ // to the left of it was pixels + remark(U[1], &upmark); + }else + upmark = assoc[U[1]]; + } + #endif + } + if(!upmark){ // there's nothing above - set current pixel to incremented counter + #ifdef LABEL_8 // check, whether pixel is not single + if( !(x && ((D && D[-1]) || ptr[-1])) // no left neighbours + && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours + && !(D && D[0])){ // no neighbour down + *ptr = 0; // erase this hermit! + continue; + } + #endif + upmark = ++Ncur; + } + + // now remark DL and DR corners (bottom pixel will be checked on next line) + if(D){ + #ifdef LABEL_8 + if(x && D[-1]){ + remark(D[-1], &upmark); + } + if(x < w && D[1]){ + remark(D[1], &upmark); + } + #else + if(D[0]) remark(D[0], &upmark); + #endif + } + // remark current + remark(curval, &upmark); + } + } +printf("time for step 2: %gs, found %zd objects\n", dtime()-t0, Ncur); +t0 = dtime(); + // Step 3: rename markers + DBG("rename markers\n"); + // first correct complex assotiations in assoc + OMP_FOR() + for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + size_t p = *ptr; + if(!p){continue;} + *ptr = assoc[p]; + } + } +printf("time for step 3: %gs\n", dtime()-t0); +#ifdef EBUG +printf("\n\n"); +for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + if(*ptr) printf("%02zx", *ptr); + else printf(" "); + } + printf("\n"); +} +#endif + FREE(assoc); + FREE(labels); diff --git a/erosion-dilation/cclabling_1.h b/erosion-dilation/cclabling_1.h new file mode 100644 index 0000000..e8a1003 --- /dev/null +++ b/erosion-dilation/cclabling_1.h @@ -0,0 +1,170 @@ +/* + * cclabling_1.h - inner part of functions cclabel4_1 and cclabel8_1 + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +double t0 = dtime(); + + CCbox *ret = NULL; + size_t N = 0, // current label + Ntot = 0, // number of found objects + Nmax = W*H/4; // max number of labels + size_t *labels = char2st(I, W, H, W_0); + int w = W - 1; +printf("time for char2st: %gs\n", dtime()-t0); + int y; + size_t *assoc = Malloc(Nmax, sizeof(size_t)); // allocate memory for "remark" array + inline void remark(size_t old, size_t *newv){ + size_t new = *newv; + if(old == new || assoc[old] == new || assoc[new] == old){ + DBG("(%zd)\n", new); + return; + } + DBG("[cur: %zx, tot: %zd], %zx -> %zx ", N, Ntot, old, new); + if(!assoc[old]){ // try to remark non-marked value + assoc[old] = new; + DBG("\n"); + return; + } + // value was remarked -> we have to remark "new" to assoc[old] + // and decrement all marks, that are greater than "new" + size_t ao = assoc[old]; + DBG("(remarked to %zx) ", ao); + DBG(" {old: %zd, new: %zd, a[o]=%zd, a[n]=%zd} ", old, new, assoc[old], assoc[new]); + // now we must check assoc[old]: if it is < new -> change *newv, else remark assoc[old] + if(ao < new) + *newv = ao; + else{ + size_t x = ao; ao = new; new = x; // swap values + } + int m = (old > new) ? old : new; + int xx = m + W / 2; + if(xx > Nmax) xx = Nmax; + DBG(" [[xx=%d]] ", xx); + OMP_FOR() + for(int i = 1; i <= xx; i++){ + size_t m = assoc[i]; + if(m < new) continue; + if(m == new){ + assoc[i] = ao; + DBG("remark %x (%zd) to %zx ", i, m, ao); + }else{ + assoc[i]--; + DBG("decr %x: %zx, ", i, assoc[i]); + } + } + DBG("\n"); + Ntot--; + } + t0 = dtime(); + size_t *ptr = labels; + for(y = 0; y < H; y++){ // FIXME!!!! throw out that fucking "if" checking coordinates!!! + bool found = false; + size_t curmark = 0; // mark of pixel to the left + for(int x = 0; x < W; x++, ptr++){ + size_t curval = *ptr; + if(!curval){ found = false; continue;} + size_t *U = (y) ? &ptr[-W] : NULL; + size_t upmark = 0; // mark from line above + if(!found){ // new pixel, check neighbours above + found = true; + // now check neighbours in upper string: + if(U){ + //#ifdef LABEL_8 + if(x && U[-1]){ // there is something in upper left corner -> use its value + upmark = U[-1]; + }else // check point above only if there's nothing in left up + //#endif + if(U[0]) upmark = U[0]; + //#ifdef LABEL_8 + if(x < w && U[1]){ // there's something in upper right + if(upmark){ // to the left of it was pixels + remark(U[1], &upmark); + }else + upmark = U[1]; + } + //#endif + } + if(!upmark){ // there's nothing above - set current pixel to incremented counter + #ifdef LABEL_8 // check, whether pixel is not single + size_t *D = (y < w) ? &ptr[W] : NULL; + if( !(x && ((D && D[-1]) || ptr[-1])) // no left neighbours + && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours + && !(D && D[0])){ // no neighbour down + *ptr = 0; // erase this hermit! + continue; + } + #endif + upmark = ++N; + Ntot++; + assoc[upmark] = upmark; // refresh "assoc" + } + *ptr = upmark; + curmark = upmark; + //remark(curval, &upmark); + }else{ // there was something to the left -> we must chek only U[1] + if(U){ + if(x < w && U[1]){ // there's something in upper right + remark(U[1], &curmark); + } + } + *ptr = curmark; + } + } + } +printf("time for step 2: %gs, found %zd objects\n", dtime()-t0, Ntot); + DBG("Initial mark\n"); + #ifdef EBUG +for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + if(*ptr) printf("%02zx", *ptr); + else printf(" "); + } + printf("\n"); +} + #endif + +t0 = dtime(); + // Step 2: rename markers + DBG("rename markers\n"); + // first correct complex assotiations in assoc + OMP_FOR() + for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + size_t p = *ptr; + if(!p){continue;} + *ptr = assoc[p]; + } + } +printf("time for step 3: %gs\n", dtime()-t0); +#ifdef EBUG +printf("\n\n"); +for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + if(*ptr) printf("%02zx", *ptr); + else printf(" "); + } + printf("\n"); +} +#endif + FREE(assoc); + FREE(labels); diff --git a/erosion-dilation/dilation.h b/erosion-dilation/dilation.h new file mode 100644 index 0000000..9489304 --- /dev/null +++ b/erosion-dilation/dilation.h @@ -0,0 +1,75 @@ +/* + * dilation.h - inner part of function `dilation` + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +/* + * HERE'S NO ANY "FILE-GUARDS" BECAUSE FILE IS MULTIPLY INCLUDED! + * I use this because don't want to dive into depths of BOOST + * + * multi-including with different defines before - is a simplest way to + * modify simple code for avoiding extra if-then-else + */ + + +#if ! defined IM_UP && ! defined IM_DOWN // image without upper and lower borders +OMP_FOR() +for(y = 1; y < h; y++) +#endif +{ + unsigned char *iptr = &image[W*y]; + unsigned char *optr = &ret[W*y]; + unsigned char p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + int x; + if(iptr[1] & 0x80) p |= 1; + *optr = p; + optr++; iptr++; + for(x = 1; x < w; x++, iptr++, optr++){ + p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + if(iptr[-1] & 1) p |= 0x80; + *optr = p; + } + p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + if(iptr[-1] & 1) p |= 0x80; + *optr = p; + optr++; iptr++; +} + diff --git a/erosion-dilation/fc_filter.h b/erosion-dilation/fc_filter.h new file mode 100644 index 0000000..504b47a --- /dev/null +++ b/erosion-dilation/fc_filter.h @@ -0,0 +1,67 @@ +/* + * fc_filter.h - inner part of function `FC_filter` + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +// The same shit as for dilation.h + +#if ! defined IM_UP && ! defined IM_DOWN +OMP_FOR() +for(y = 1; y < h; y++) +#endif +{ + unsigned char *iptr = &image[W*y]; + unsigned char *optr = &ret[W*y]; + unsigned char p = *iptr << 1 | *iptr >> 1 + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + int x; + if(iptr[1] & 0x80) p |= 1; + *optr = p & *iptr; + optr++; iptr++; + for(x = 1; x < w; x++, iptr++, optr++){ + p = *iptr << 1 | *iptr >> 1 + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + if(iptr[-1] & 1) p |= 0x80; + *optr = p & *iptr; + } + p = *iptr << 1 | *iptr >> 1 + #ifndef IM_UP + | iptr[-W] + #endif + #ifndef IM_DOWN + | iptr[W] + #endif + ; + if(iptr[-1] & 1) p |= 0x80; + *optr = p & *iptr; + optr++; iptr++; +} diff --git a/erosion-dilation/fs_filter.h b/erosion-dilation/fs_filter.h new file mode 100644 index 0000000..e69de29 diff --git a/erosion-dilation/main.c b/erosion-dilation/main.c new file mode 100644 index 0000000..3aaff7a --- /dev/null +++ b/erosion-dilation/main.c @@ -0,0 +1,173 @@ +/* + * main.c - test file for binmorph.c + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include "binmorph.h" + +// these includes are for randini +#include +#include +#include +#include +/* + * Generate a quasy-random number to initialize PRNG + * name: throw_random_seed + * @return value for curandSetPseudoRandomGeneratorSeed or srand48 + */ +long throw_random_seed(){ + //FNAME(); + long r_ini; + int fail = 0; + int fd = open("/dev/random", O_RDONLY); + do{ + if(-1 == fd){ + fail = 1; break; + } + if(sizeof(long) != read(fd, &r_ini, sizeof(long))){ + fail = 1; + } + close(fd); + }while(0); + if(fail){ + double tt = dtime() * 1e6; + double mx = (double)LONG_MAX; + r_ini = (long)(tt - mx * floor(tt/mx)); + } + return (r_ini); +} + + +int main(int ac, char **v){ + int i,j, W = 28, H = 28, W_0; + //double midX = (W - 1.0)/ 4. - 1., midY = (H - 1.) / 4. - 1.; + //double ro = sqrt(midX*midY), ri = ro / 1.5; + bool *inp = Malloc(W * H, sizeof(bool)); + printf("\n"); + srand48(throw_random_seed()); + for(i = 0; i < H; i++) + for(j = 0; j < W; j++) + inp[i*W+j] = (drand48() > 0.5); + /* + for(i = 0; i < H/2; i++){ + for(j = 0; j < W/2; j++){ + double x = j - midX, y = i - midY; + double r = sqrt(x*x + y*y); + if((r < ro && r > ri) || fabs(x) == fabs(y)) + inp[i*W+j] = inp[i*W+j+W/2] = inp[(i+H/2)*W+j] = inp[(i+H/2)*W+j+W/2] = true; + } + }*/ + unsigned char *b2ch = bool2char(inp, W, H, &W_0); + FREE(inp); + printf("inpt: (%dx%d)\n", W_0, H); + printC(b2ch, W_0, H); + unsigned char *eros = erosion(b2ch, W_0, H); + printf("erosion: (%dx%d)\n", W_0, H); + printC(eros, W_0, H); + unsigned char *dilat = dilation(b2ch, W_0, H); + printf("dilation: (%dx%d)\n", W_0, H); + printC(dilat, W_0, H); + unsigned char *erdilat = dilation(eros, W_0, H); + printf("\n\ndilation of erosion (openinfg: (%dx%d)\n", W_0, H); + printC(erdilat, W_0, H); + unsigned char *dileros = erosion(dilat, W_0, H); + printf("erosion of dilation (closing): (%dx%d)\n", W_0, H); + printC(dileros, W_0, H); + printf("\n\n\n image - opening of original minus original:\n"); + unsigned char *immer = substim(b2ch, erdilat, W_0, H); + printC(immer, W_0, H); + FREE(eros); FREE(dilat); FREE(erdilat); FREE(dileros); + inp = char2bool(immer, W, H, W_0); + printf("\n\nAnd boolean for previous image (%dx%d):\n ",W,H); + printB(inp, W, H); + FREE(immer); + immer = FC_filter(b2ch, W_0, H); + printf("\n\n\nFilter for 4-connected areas searching:\n"); + printC(immer, W_0, H); + FREE(immer); + size_t NL; + printf("\nmark 8-connected components:\n"); + cclabel8_1(b2ch, W, H, W_0, &NL); + printf("\nmark 8-connected components (old):\n"); + cclabel8(b2ch, W, H, W_0, &NL); + FREE(b2ch); + + return 0; + W = H = 1000; + FREE(inp); + inp = Malloc(W * H, sizeof(bool)); + for(i = 0; i < H; i++) + for(j = 0; j < W; j++) + inp[i*W+j] = (drand48() > 0.5); + b2ch = bool2char(inp, W, H, &W_0); + FREE(inp); + printf("\n\n\n1) 1 cc4 for 2000x2000 .. "); + fflush(stdout); + double t0 = dtime(); + for(i = 0; i < 1; i++){ + CCbox *b = cclabel4(b2ch, W, H, W_0, &NL); + FREE(b); + } + printf("%.3f s\n", dtime() - t0); + printf("\n2) 1 cc8 for 2000x2000 .. "); + fflush(stdout); + t0 = dtime(); + for(i = 0; i < 1; i++){ + CCbox *b = cclabel8(b2ch, W, H, W_0, &NL); + FREE(b); + } + printf("%.3f s\n", dtime() - t0); + printf("\n3) 1 cc8_1 for 2000x2000 .. "); + fflush(stdout); + t0 = dtime(); + for(i = 0; i < 1; i++){ + CCbox *b = cclabel8_1(b2ch, W, H, W_0, &NL); + FREE(b); + } + printf("%.3f s\n", dtime() - t0); + + /* printf("\n\n\n\nSome tests:\n1) 10000 dilations for 2000x2000 .. "); + fflush(stdout); + double t0 = dtime(); + for(i = 0; i < 10000; i++){ + unsigned char *dilat = dilation(b2ch, W, H); + free(dilat); + } + printf("%.3f s\n", dtime() - t0); + printf("2) 10000 erosions for 2000x2000 .. "); + fflush(stdout); + t0 = dtime(); + for(i = 0; i < 10000; i++){ + unsigned char *dilat = erosion(b2ch, W, H); + free(dilat); + } + printf("%.3f s\n", dtime() - t0); + printf("3) 10000 image substitutions for 2000x2000 .. "); + fflush(stdout); + unsigned char *dilat = dilation(b2ch, W, H); + t0 = dtime(); + for(i = 0; i < 10000; i++){ + unsigned char *res = substim(b2ch, dilat, W, H); + free(res); + } + printf("%.3f s\n", dtime() - t0);*/ + return 0; +} diff --git a/fifo_lifo.c b/fifo_lifo.c new file mode 100644 index 0000000..b429c02 --- /dev/null +++ b/fifo_lifo.c @@ -0,0 +1,129 @@ +/* + * fifo_lifo.c - simple FIFO/LIFO + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include +#include + +#include "fifo_lifo.h" + +#ifndef Malloc +#define Malloc my_alloc +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} +#endif +#ifndef FREE +#define FREE(arg) do{free(arg); arg = NULL;}while(0) +#endif + + +/** + * push data into the tail of a stack (like FIFO) + * @param lst (io) - list + * @param v (i) - data to push + * @return pointer to just pused node + */ +List *push_tail(List **lst, Ldata v){ + List *node; + if(!lst) return NULL; + if((node = (List*) Malloc(1, sizeof(List))) == 0) + return NULL; // allocation error + node->data = v; // insert data + if(!*lst){ + *lst = node; + (*lst)->last = node; + return node; + } + (*lst)->last->next = node; + (*lst)->last = node; + return node; +} + +/** + * push data into the head of a stack (like LIFO) + * @param lst (io) - list + * @param v (i) - data to push + * @return pointer to just pused node + */ +List *push(List **lst, Ldata v){ + List *node; + if(!lst) return NULL; + if((node = (List*) Malloc(1, sizeof(List))) == 0) + return NULL; // allocation error + node->data = v; // insert data + if(!*lst){ + *lst = node; + (*lst)->last = node; + return node; + } + node->next = *lst; + node->last = (*lst)->last; + *lst = node; + return node; +} + +/** + * get data from head of list + * @param lst (io) - list + * @return data from lst head + */ +Ldata pop(List **lst){ + Ldata ret; + List *node = *lst; + if(!lst) errx(1, "NULL pointer to buffer"); + if(!*lst) errx(1, "Empty buffer"); + ret = (*lst)->data; + *lst = (*lst)->next; + FREE(node); + return ret; +} + +#ifdef STANDALONE +int main(void) { + List *f = NULL; + int i, d, ins[] = {4,2,6,1,3,4,7}; + for(i = 0; i < 7; i++) + if(!(push(&f, ins[i]))) + err(1, "Malloc error"); // can't insert + else printf("pushed %d\n", ins[i]); + while(f){ + d = pop(&f); + printf("pulled: %d\n", d); + } + for(i = 0; i < 7; i++) + if(!(push_tail(&f, ins[i]))) + err(1, "Malloc error"); // can't insert + else printf("pushed to head %d\n", ins[i]); + while(f){ + d = pop(&f); + printf("pulled: %d\n", d); + } + return 0; +} +#endif diff --git a/fifo_lifo.h b/fifo_lifo.h new file mode 100644 index 0000000..0eb79cc --- /dev/null +++ b/fifo_lifo.h @@ -0,0 +1,37 @@ +/* + * fifo_lifo.h + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __FIFO_LIFO_H__ +#define __FIFO_LIFO_H__ + +typedef int Ldata; + +typedef struct buff_node{ + Ldata data; + struct buff_node *next, *last; +} List; + +List *push_tail(List **lst, Ldata v); +List *push(List **lst, Ldata v); +Ldata pop(List **lst); + +#endif // __FIFO_LIFO_H__ diff --git a/getopt/README b/getopt/README new file mode 100644 index 0000000..5b40722 --- /dev/null +++ b/getopt/README @@ -0,0 +1,9 @@ +1. Functions for simplify control of complex parameters +cmdlnopts.c - example of use +parceargs.c - main functions +parceargs.h - their header + +2. Functions for regular getopt_long +getopt.c - snippet file +getopt.h - its header + diff --git a/getopt/cmdlnopts.c b/getopt/cmdlnopts.c new file mode 100644 index 0000000..8577e08 --- /dev/null +++ b/getopt/cmdlnopts.c @@ -0,0 +1,123 @@ +/* + * cmdlnopts.c - the only function that parce cmdln args and returns glob parameters + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include "cmdlnopts.h" + +/* + * here are global parameters initialisation + */ + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +myoption cmdlnopts[] = { + {"help", 0, NULL, 'h', arg_int, APTR(&help), N_("show this help")}, + // ... + end_option +}; + + +/** + * Parse string of suboptions (--option=name1=var1:name2=var2... or -O name1=var1,name2=var2...) + * Suboptions could be divided by colon or comma + * + * !!NAMES OF SUBOPTIONS ARE CASE-UNSENSITIVE!!! + * + * @param arg (i) - string with description + * @param V (io) - pointer to suboptions array (result will be stored in sopts->val) + * @return TRUE if success + */ +bool get_suboptions(void *arg, suboptions *V){ + char *tok, *val, *par; + int i; + tok = strtok(arg, ":,"); + do{ + if((val = strchr(tok, '=')) == NULL){ // wrong format + WARNX(_("Wrong format: no value for keyword")); + return FALSE; + } + *val++ = '\0'; + par = tok; + for(i = 0; V[i].val; i++){ + if(strcasecmp(par, V[i].par) == 0){ // found parameter + if(V[i].isdegr){ // DMS + if(!get_radians(V[i].val, val)) // wrong angle + return FALSE; + DBG("Angle: %g rad\n", *(V[i].val)); + }else{ // simple float + if(!myatof(V[i].val, val)) // wrong number + return FALSE; + DBG("Float val: %g\n", *(V[i].val)); + } + break; + } + } + if(!V[i].val){ // nothing found - wrong format + WARNX(_("Bad keyword!")); + return FALSE; + } + }while((tok = strtok(NULL, ":,"))); + return TRUE; +} + +/** + * functions of subargs parcing can looks as this + */ +bool get_mir_par(void *arg, int N _U_){ + suboptions V[] = { // array of mirror parameters and string keys for cmdln pars + {&M.D, "diam", FALSE}, + {&M.F, "foc", FALSE}, + {&M.Zincl, "zincl", TRUE}, + {&M.Aincl, "aincl", TRUE}, + {&M.objA, "aobj", TRUE}, + {&M.objZ, "zobj", TRUE}, + {&M.foc, "ccd", FALSE}, + {0,0,0} + }; + return get_suboptions(arg, V); +} + +/** + * Parce command line options and return dynamically allocated structure + * to global parameters + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +glob_pars *parce_args(int argc, char **argv){ + int i; + void *ptr; + ptr = memcpy(&G, &Gdefault, sizeof(G)); assert(ptr); + ptr = memcpy(&M, &Mdefault, sizeof(M)); assert(ptr); + G.Mirror = &M; + // format of help: "Usage: progname [args]\n" + change_helpstring("Usage: %s [args]\n\n\tWhere args are:\n"); + // parse arguments + parceargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(argc > 0){ + printf("\nIgnore argument[s]:\n"); + for (i = 0; i < argc; i++) + printf("\t%s\n", argv[i]); + } + return &G; +} + diff --git a/getopt/cmdlnopts.h b/getopt/cmdlnopts.h new file mode 100644 index 0000000..4778788 --- /dev/null +++ b/getopt/cmdlnopts.h @@ -0,0 +1,36 @@ +/* + * cmdlnopts.h - comand line options for parceargs + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __CMDLNOPTS_H__ +#define __CMDLNOPTS_H__ + +#include "parceargs.h" + +/* + * here are some typedef's for global data + */ + +typedef SomeType glob_pars; + +glob_pars *parce_args(int argc, char **argv); + +#endif // __CMDLNOPTS_H__ diff --git a/getopt/getopt.c b/getopt/getopt.c new file mode 100644 index 0000000..c4dd475 --- /dev/null +++ b/getopt/getopt.c @@ -0,0 +1,128 @@ +/* + * getopt.c - simple functions for getopt_long + * provide functions + * usage - to show help message & exit + * parse_args - to parce argv & fill global flags + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include +#include + + +/* + * here are global variables and global data structures initialisation, like + * int val = 10; // default value + */ + +/* + * HERE A PART of main.c: + * + setlocale(LC_ALL, ""); + setlocale(LC_NUMERIC, "C"); + bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR); + textdomain(GETTEXT_PACKAGE); + parse_args(argc, argv); // <- set all flags + * + * ITS HEADER: +#ifndef GETTEXT_PACKAGE // should be defined by make +#define GETTEXT_PACKAGE "progname" +#endif +#ifndef LOCALEDIR // should be defined by make +#define LOCALEDIR "/path/to/locale" +#endif + */ + +void usage(char *fmt, ...){ + va_list ap; + va_start(ap, fmt); + printf("\n"); + if (fmt != NULL){ + vprintf(fmt, ap); + printf("\n\n"); + } + va_end(ap); + // ":\t%s [] [ ]\n" + printf(_("Usage:\t%s [options] [output files prefix]\n"), + __progname); + // "\t:\n" + printf(_("\tOptions:\n")); + printf("\t-A,\t--author=author\t\t%s\n", + // " " + _("program author")); + // and so on + exit(1); +} + +void parse_args(int argc, char **argv){ + int i; + char short_options[] = "A"; // all short equivalents + struct option long_options[] = { +/* { name, has_arg, flag, val }, where: + * name - name of long parameter + * has_arg = 0 - no argument, 1 - need argument, 2 - unnesessary argument + * flag = NULL to return val, pointer to int - to set it + * value of val (function returns 0) + * val - getopt_long return value or value, flag setting to + * !!! last string - for zeros !!! + */ + {"author", 1, 0, 'A'}, + // and so on + { 0, 0, 0, 0 } + }; + if(argc == 1){ + // " " + usage(_("Any parameters are absent")); + } + while(1){ + int opt; + if((opt = getopt_long(argc, argv, short_options, + long_options, NULL)) == -1) break; + switch(opt){ + case 0: // only long option + // do something? + break; + case 'A': + author = strdup(optarg); + // " : %s" + //info(_("Program author: %s"), author); + break; + // ... + default: + usage(NULL); + } + } + argc -= optind; + argv += optind; + if(argc == 0){ + // there's no free parameters + }else{ + /* there was free parameter + * for example, filename + outfile = argv[0]; + argc--; + argv++; + */ + } + if(argc > 0){ + // " []:\n" + printf(_("Ignore argument[s]:\n")); + for (i = 0; i < argc; i++) + warnx("%s ", argv[i]); + } +} diff --git a/getopt/getopt.h b/getopt/getopt.h new file mode 100644 index 0000000..d913e3f --- /dev/null +++ b/getopt/getopt.h @@ -0,0 +1,18 @@ +#ifndef __USAGE_H__ +#define __USAGE_H__ + +#include "takepic.h" +#include +#include + +/* + * here are global variables and global data structures definition, like + * int val; + * enum{reset = 0, set}; + */ + +extern char *__progname; +void usage(char *fmt, ...); +void parse_args(int argc, char **argv); + +#endif // __USAGE_H__ diff --git a/getopt/parceargs.c b/getopt/parceargs.c new file mode 100644 index 0000000..d7b84ca --- /dev/null +++ b/getopt/parceargs.c @@ -0,0 +1,316 @@ +/* + * parceargs.c - parcing command line arguments & print help + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include // DBG +#include // getopt_long +#include // calloc, exit, strtoll +#include // assert +#include // strdup, strchr, strlen +#include // INT_MAX & so on +#include // gettext +#include // isalpha +#include "parceargs.h" + +#ifdef EBUG + #define DBG(...) printf(__VA_ARGS__) +#else + #define DBG(...) +#endif + +// macro to print help messages +#ifndef PRNT + #define PRNT(x) gettext(x) +#endif + +char *helpstring = "%s\n"; + +/** + * Change standard help header + * MAY consist ONE "%s" for progname + * @param str (i) - new format + */ +void change_helpstring(char *s){ + int pcount = 0, scount = 0; + char *str = s; + // check `helpstring` and set it to default in case of error + for(; pcount < 2; str += 2){ + if(!(str = strchr(str, '%'))) break; + if(str[1] != '%') pcount++; // increment '%' counter if it isn't "%%" + else{ + str += 2; // pass next '%' + continue; + } + if(str[1] == 's') scount++; // increment "%s" counter + }; + DBG("pc: %d, sc: %d\n", pcount, scount); + if(pcount > 1 || pcount != scount){ // amount of pcount and/or scount wrong + fprintf(stderr, "Wrong helpstring!\n"); + exit(-1); + } + helpstring = s; + DBG("hs: %s\n", helpstring); +} + +/** + * Carefull atoll/atoi + * @param num (o) - returning value (or NULL if you wish only check number) - allocated by user + * @param str (i) - string with number must not be NULL + * @param t (i) - T_INT for integer or T_LLONG for long long (if argtype would be wided, may add more) + * @return TRUE if conversion sone without errors, FALSE otherwise + */ +bool myatoll(void *num, char *str, argtype t){ + long long tmp, *llptr; + int *iptr; + char *endptr; + assert(str); + assert(num); + tmp = strtoll(str, &endptr, 0); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_longlong: + llptr = (long long*) num; + *llptr = tmp; + break; + case arg_int: + default: + if(tmp < INT_MIN || tmp > INT_MAX){ + fprintf(stderr, "Integer out of range\n"); + return FALSE; + } + iptr = (int*)num; + *iptr = (int)tmp; + } + return TRUE; +} + +// the same as myatoll but for double +// There's no NAN & INF checking here (what if they would be needed?) +bool myatod(void *num, const char *str, argtype t){ + double tmp, *dptr; + float *fptr; + char *endptr; + assert(str); + tmp = strtod(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_double: + dptr = (double *) num; + *dptr = tmp; + break; + case arg_float: + default: + fptr = (float *) num; + *fptr = (float)tmp; + break; + } + return TRUE; +} + +/** + * Get index of current option in array options + * @param opt (i) - returning val of getopt_long + * @param options (i) - array of options + * @return index in array + */ +int get_optind(int opt, myoption *options){ + int oind; + myoption *opts = options; + assert(opts); + for(oind = 0; opts->name && opts->val != opt; oind++, opts++); + if(!opts->name || opts->val != opt) // no such parameter + showhelp(-1, options); + return oind; +} + +/** + * Parce command line arguments + * ! If arg is string, then value will be strdup'ed! + * + * @param argc (io) - address of argc of main(), return value of argc stay after `getopt` + * @param argv (io) - address of argv of main(), return pointer to argv stay after `getopt` + * BE CAREFUL! if you wanna use full argc & argv, save their original values before + * calling this function + * @param options (i) - array of `myoption` for arguments parcing + * + * @exit: in case of error this function show help & make `exit(-1)` + */ +void parceargs(int *argc, char ***argv, myoption *options){ + char *short_options, *soptr; + struct option *long_options, *loptr; + size_t optsize, i; + myoption *opts = options; + // check whether there is at least one options + assert(opts); + assert(opts[0].name); + // first we count how much values are in opts + for(optsize = 0; opts->name; optsize++, opts++); + // now we can allocate memory + short_options = calloc(optsize * 3 + 1, 1); // multiply by three for '::' in case of args in opts + long_options = calloc(optsize + 1, sizeof(struct option)); + opts = options; loptr = long_options; soptr = short_options; + // fill short/long parameters and make a simple checking + for(i = 0; i < optsize; i++, loptr++, opts++){ + // check + assert(opts->name); // check name + if(opts->has_arg){ + assert(opts->type != arg_none); // check error with arg type + assert(opts->argptr); // check pointer + } + if(opts->type != arg_none) // if there is a flag without arg, check its pointer + assert(opts->argptr); + // fill long_options + // don't do memcmp: what if there would be different alignment? + loptr->name = opts->name; + loptr->has_arg = opts->has_arg; + loptr->flag = opts->flag; + loptr->val = opts->val; + // fill short options if they are: + if(!opts->flag){ + *soptr++ = opts->val; + if(opts->has_arg) // add ':' if option has required argument + *soptr++ = ':'; + if(opts->has_arg == 2) // add '::' if option has optional argument + *soptr++ = ':'; + } + } + // now we have both long_options & short_options and can parse `getopt_long` + while(1){ + int opt; + int oindex = 0, optind = 0; // oindex - number of option in argv, optind - number in options[] + if((opt = getopt_long(*argc, *argv, short_options, long_options, &oindex)) == -1) break; + if(opt == '?'){ + opt = optopt; + optind = get_optind(opt, options); + if(options[optind].has_arg == 1) showhelp(optind, options); // need argument + } + else{ + if(opt == 0 || oindex > 0) optind = oindex; + else optind = get_optind(opt, options); + } + opts = &options[optind]; +DBG ("\n*******\noption %s (oindex = %d / optind = %d)", options[optind].name, oindex, optind); +if(optarg) DBG (" with arg %s", optarg); +DBG ("\n"); + if(opt == 0 && opts->has_arg == 0) continue; // only long option changing integer flag +DBG("opt = %c, arg type: ", opt); + // now check option + if(opts->has_arg == 1) assert(optarg); + bool result = TRUE; + // even if there is no argument, but argptr != NULL, think that optarg = "1" + if(!optarg) optarg = "1"; + switch(opts->type){ + default: + case arg_none: +DBG("none\n"); + if(opts->argptr) *((int*)opts->argptr) = 1; // set argptr to 1 + break; + case arg_int: +DBG("integer\n"); + result = myatoll(opts->argptr, optarg, arg_int); + break; + case arg_longlong: +DBG("long long\n"); + result = myatoll(opts->argptr, optarg, arg_longlong); + break; + case arg_double: +DBG("double\n"); + result = myatod(opts->argptr, optarg, arg_double); + break; + case arg_float: +DBG("double\n"); + result = myatod(opts->argptr, optarg, arg_float); + break; + case arg_string: +DBG("string\n"); + result = (*((char **)opts->argptr) = strdup(optarg)); + break; + case arg_function: +DBG("function\n"); + result = ((argfn)opts->argptr)(optarg, optind); + break; + } + if(!result){ +DBG("OOOPS! Error in result\n"); + showhelp(optind, options); + } + } + *argc -= optind; + *argv += optind; +} + +/** + * Show help information based on myoption->help values + * @param oindex (i) - if non-negative, show only help by myoption[oindex].help + * @param options (i) - array of `myoption` + * + * @exit: run `exit(-1)` !!! + */ +void showhelp(int oindex, myoption *options){ + // ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext, + // but you can redefine it before `#include "parceargs.h"` + int max_opt_len = 0; // max len of options substring - for right indentation + const int bufsz = 255; + char buf[bufsz+1]; + myoption *opts = options; + assert(opts); + assert(opts[0].name); // check whether there is at least one options + if(oindex > -1){ // print only one message + opts = &options[oindex]; + printf(" "); + if(!opts->flag && isalpha(opts->val)) printf("-%c, ", opts->val); + printf("--%s", opts->name); + if(opts->has_arg == 1) printf("=arg"); + else if(opts->has_arg == 2) printf("[=arg]"); + printf(" %s\n", PRNT(opts->help)); + exit(-1); + } + // header, by default is just "progname\n" + printf("\n"); + if(strstr(helpstring, "%s")) // print progname + printf(helpstring, __progname); + else // only text + printf("%s", helpstring); + printf("\n"); + // count max_opt_len + do{ + int L = strlen(opts->name); + if(max_opt_len < L) max_opt_len = L; + }while((++opts)->name); + max_opt_len += 14; // format: '-S , --long[=arg]' - get addition 13 symbols + opts = options; + // Now print all help + do{ + int p = sprintf(buf, " "); // a little indent + if(!opts->flag && isalpha(opts->val)) // .val is short argument + p += snprintf(buf+p, bufsz-p, "-%c, ", opts->val); + p += snprintf(buf+p, bufsz-p, "--%s", opts->name); + if(opts->has_arg == 1) // required argument + p += snprintf(buf+p, bufsz-p, "=arg"); + else if(opts->has_arg == 2) // optional argument + p += snprintf(buf+p, bufsz-p, "[=arg]"); + assert(p < max_opt_len); // there would be magic if p >= max_opt_len + printf("%-*s%s\n", max_opt_len+1, buf, PRNT(opts->help)); // write options & at least 2 spaces after + }while((++opts)->name); + printf("\n\n"); + exit(-1); +} diff --git a/getopt/parceargs.h b/getopt/parceargs.h new file mode 100644 index 0000000..554803b --- /dev/null +++ b/getopt/parceargs.h @@ -0,0 +1,105 @@ +/* + * parceargs.h - headers for parcing command line arguments + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __PARCEARGS_H__ +#define __PARCEARGS_H__ + +#include // bool +#include + +#ifndef TRUE + #define TRUE true +#endif + +#ifndef FALSE + #define FALSE false +#endif + +// macro for argptr +#define APTR(x) ((void*)x) + +// if argptr is a function: +typedef bool(*argfn)(void *arg, int N); + +/* + * type of getopt's argument + * WARNING! + * My function change value of flags by pointer, so if you want to use another type + * make a latter conversion, example: + * char charg; + * int iarg; + * myoption opts[] = { + * {"value", 1, NULL, 'v', arg_int, &iarg, "char val"}, ..., end_option}; + * ..(parce args).. + * charg = (char) iarg; + */ +typedef enum { + arg_none = 0, // no arg + arg_int, // integer + arg_longlong, // long long + arg_double, // double + arg_float, // float + arg_string, // char * + arg_function // parce_args will run function `bool (*fn)(char *optarg, int N)` +} argtype; + +/* + * Structure for getopt_long & help + * BE CAREFUL: .argptr is pointer to data or pointer to function, + * conversion depends on .type + * + * ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext, + * but you can redefine it before `#include "parceargs.h"` + * + * if arg is string, then value wil be strdup'ed like that: + * char *str; + * myoption opts[] = {{"string", 1, NULL, 's', arg_string, &str, "string val"}, ..., end_option}; + * *(opts[1].str) = strdup(optarg); + * in other cases argptr should be address of some variable (or pointer to allocated memory) + * + * NON-NULL argptr should be written inside macro APTR(argptr) or directly: (void*)argptr + * + * !!!LAST VALUE OF ARRAY SHOULD BE `end_option` or ZEROS !!! + * + */ +typedef struct{ + // these are from struct option: + const char *name; // long option's name + int has_arg; // 0 - no args, 1 - nesessary arg, 2 - optionally arg + int *flag; // NULL to return val, pointer to int - to set its value of val (function returns 0) + int val; // short opt name (if flag == NULL) or flag's value + // and these are mine: + argtype type; // type of argument + void *argptr; // pointer to variable to assign optarg value or function `bool (*fn)(char *optarg, int N)` + char *help; // help string which would be shown in function `showhelp` or NULL +} myoption; + +// last string of array (all zeros) +#define end_option {0,0,0,0,0,0,0} + + +extern const char *__progname; + +void showhelp(int oindex, myoption *options); +void parceargs(int *argc, char ***argv, myoption *options); +void change_helpstring(char *s); + +#endif // __PARCEARGS_H__ diff --git a/kicad-copy/Makefile b/kicad-copy/Makefile new file mode 100644 index 0000000..33c1573 --- /dev/null +++ b/kicad-copy/Makefile @@ -0,0 +1,23 @@ +PROGRAM = copy_pcb +LDFLAGS = +SRCS = cmdlnopts.c copy_pcb.c parceargs.c usefull_macros.c +CC = gcc +DEFINES = -D_XOPEN_SOURCE=1001 -DGETTEXT_PACKAGE=\"copy_pcb\" -DLOCALEDIR=\".\" +DEFINES +=-DEBUG +CXX = gcc +CFLAGS = -std=gnu99 -Wall -Werror $(DEFINES) +OBJS = $(SRCS:.c=.o) +all : $(PROGRAM) +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +# some addition dependencies +# %.o: %.c +# $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@ +#$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS) +# @touch $@ + +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(SRCS) diff --git a/kicad-copy/cmdlnopts.c b/kicad-copy/cmdlnopts.c new file mode 100644 index 0000000..676730a --- /dev/null +++ b/kicad-copy/cmdlnopts.c @@ -0,0 +1,65 @@ +/* + * cmdlnopts.c - the only function that parce cmdln args and returns glob parameters + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#include "copy_pcb.h" +#include "cmdlnopts.h" + +/* + * here are global parameters initialisation + */ +glob_pars G; +int help = 0; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +myoption cmdlnopts[] = { + {"help", 0, NULL, 'h', arg_int, APTR(&help), N_("show this help")}, + {"format", 1, NULL, 'f', arg_string, APTR(&G.format), N_("printf-like format of similar items (%s - MUST BE - old name; %n - new num)")}, + {"nparts", 1, NULL, 'N', arg_int, APTR(&G.nparts), N_("Number of copies")}, + {"list", 1, NULL, 'L', arg_string, APTR(&G.list), N_("comma-separated list of sheets")}, + {"input", 1, NULL, 'i', arg_string, APTR(&G.ifile), N_("input file name")}, + end_option +}; + +/** + * Parce command line options and return dynamically allocated structure + * to global parameters + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +glob_pars *parce_args(int argc, char **argv){ + int i; + memset(&G, sizeof(glob_pars), 0); // clear all + // format of help: "Usage: progname [args]\n" + change_helpstring("Usage: %s [args]\n\n\tWhere args are:\n"); + // parse arguments + parceargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(argc > 0){ + printf("\nIgnore argument[s]:\n"); + for (i = 0; i < argc; i++) + printf("\t%s\n", argv[i]); + } + return &G; +} + diff --git a/kicad-copy/cmdlnopts.h b/kicad-copy/cmdlnopts.h new file mode 100644 index 0000000..7c67f4a --- /dev/null +++ b/kicad-copy/cmdlnopts.h @@ -0,0 +1,43 @@ +/* + * cmdlnopts.h - comand line options for parceargs + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __CMDLNOPTS_H__ +#define __CMDLNOPTS_H__ + +#include "parceargs.h" + +/* + * here are some typedef's for global data + */ + +typedef struct{ + char *format; // printf-like format: "%s_%d" + int nparts; // number of parts from list + char *list; // comma-separated list + char *ifile; // input file +}glob_pars; + +extern glob_pars G; + +glob_pars *parce_args(int argc, char **argv); + +#endif // __CMDLNOPTS_H__ diff --git a/kicad-copy/copy_pcb.c b/kicad-copy/copy_pcb.c new file mode 100644 index 0000000..9135bfc --- /dev/null +++ b/kicad-copy/copy_pcb.c @@ -0,0 +1,115 @@ +/* + * copy_pcb.c - copy kicad PCBs + * + * Copyright 2014 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include "copy_pcb.h" + +char **List; + +/** + * Printout string from ptr0 (including) to ptr1 (excluding) and return ptr1 + */ +char *printdata(char *ptr0, char *ptr1){ + char *ptr = ptr0; + while(ptr < ptr1) printf("%c", *ptr++); + return ptr1; +} + +int build_list(){ + char *tok; + int L, i; + tok = G.list; + for(L = 0; strchr(tok, ','); L++, tok++); + DBG("found: %d commas\n", L); + if(G.nparts > 0 && G.nparts < L) L = G.nparts; + List = MALLOC(char*, L); + tok = strtok(G.list, ","); + if(!tok) return 0; + i = 0; + do{ + List[i++] = strdup(tok); + }while((tok = strtok(NULL, ",")) && i < L); + return i; +} + +int main(int argc, char **argv){ + int i, j, L; + mmapbuf *input; + initial_setup(); + parce_args(argc, argv); + if(!G.ifile){ + ERRX("You shold give input file name\n"); + return(-1); + } + if(G.format) DBG("Format: %s\n", G.format); + else{ + DBG("No format defined, use \"%%s.%%d\"\n"); + G.format = "%s.%d"; + } + if(G.nparts) DBG("Use first %d copies of list\n", G.nparts); + else DBG("Use all list given\n"); + if(!G.list){ + WARNX("Error! You should say\n\t\t%s -L list\n", __progname); + ERR("\twhere \"list\" is a comma-separated list of sheets\n"); + return(-1); + } + L = build_list(); + DBG("Found %d sheets in list:\n", L); + if(L == 0){ + ERRX("Error: empty list\n"); + return(-1); + } + for(i = 0; i < L; i++) DBG("\titem %d: %s\n", i, List[i]); + if(!(input = My_mmap(G.ifile))){ + ERRX("A strange error: can't mmap input file\n"); + return(-1); + } + char *iptr = input->data, *oldptr = iptr; + + i = 0; + char comp[32], val[32]; + int N1, N2; + while((iptr = strstr(iptr, "$Comp\n"))){ + DBG("%dth component: ", ++i); + iptr = strstr(iptr, "\nL "); + if(!iptr) break; iptr += 3; + if(sscanf(iptr, "%s %s\n", comp, val) != 2) continue; + DBG("component %s with label %s\n", comp, val); + iptr = strstr(iptr, "\nU "); + if(!iptr) break; iptr += 3; + if(sscanf(iptr, "%d %d %s\n",&N1,&N2,comp) != 3) continue; + DBG("N1 = %d; N2 = %d; comp label: %s\n",N1,N2,comp); + iptr = strstr(iptr, "\nF"); // go to line "F 0" + if(!iptr) break; iptr++; + // printout all what was before: + oldptr = printdata(oldptr, iptr); + for(j = 0; j < L; j++){ + printf("AR Path=\"/%s/%s\" Ref=\"", List[j], comp); + printf(G.format, val, j+1); + printf("\" Part=\"1\"\n"); + } + } + // printout the rest of file + if(!iptr) iptr = input->data + input->len; + printdata(oldptr, iptr); + My_munmap(input); + return 0; +} + diff --git a/kicad-copy/copy_pcb.h b/kicad-copy/copy_pcb.h new file mode 100644 index 0000000..3840c8e --- /dev/null +++ b/kicad-copy/copy_pcb.h @@ -0,0 +1,31 @@ +/* + * copy_pcb.h + * + * Copyright 2014 Edward V. Emelianov + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __COPY_PCB_H__ +#define __COPY_PCB_H__ + +#include "usefull_macros.h" +#include "cmdlnopts.h" + +#endif // __COPY_PCB_H__ + + diff --git a/kicad-copy/parceargs.c b/kicad-copy/parceargs.c new file mode 100644 index 0000000..2b53ac1 --- /dev/null +++ b/kicad-copy/parceargs.c @@ -0,0 +1,315 @@ +/* + * parceargs.c - parcing command line arguments & print help + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include // DBG +#include // getopt_long +#include // calloc, exit, strtoll +#include // assert +#include // strdup, strchr, strlen +#include // INT_MAX & so on +#include // gettext +#include // isalpha +#include "parceargs.h" + +#ifdef EBUG + #define DBG(...) printf(__VA_ARGS__) +#else + #define DBG(...) +#endif + +// macro to print help messages +#ifndef PRNT + #define PRNT(x) gettext(x) +#endif + +char *helpstring = "%s\n"; + +/** + * Change standard help header + * MAY consist ONE "%s" for progname + * @param str (i) - new format + */ +void change_helpstring(char *s){ + int pcount = 0, scount = 0; + char *str = s; + // check `helpstring` and set it to default in case of error + for(; pcount < 2; str += 2){ + if(!(str = strchr(str, '%'))) break; + if(str[1] != '%') pcount++; // increment '%' counter if it isn't "%%" + else{ + str += 2; // pass next '%' + continue; + } + if(str[1] == 's') scount++; // increment "%s" counter + }; + DBG("pc: %d, sc: %d\n", pcount, scount); + if(pcount > 1 || pcount != scount){ // amount of pcount and/or scount wrong + fprintf(stderr, "Wrong helpstring!\n"); + exit(-1); + } + helpstring = s; + DBG("hs: %s\n", helpstring); +} + +/** + * Carefull atoll/atoi + * @param num (o) - returning value (or NULL if you wish only check number) - allocated by user + * @param str (i) - string with number must not be NULL + * @param t (i) - T_INT for integer or T_LLONG for long long (if argtype would be wided, may add more) + * @return TRUE if conversion sone without errors, FALSE otherwise + */ +bool myatoll(void *num, char *str, argtype t){ + long long tmp, *llptr; + int *iptr; + char *endptr; + assert(str); + assert(num); + tmp = strtoll(str, &endptr, 0); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_longlong: + llptr = (long long*) num; + *llptr = tmp; + break; + case arg_int: + default: + if(tmp < INT_MIN || tmp > INT_MAX){ + fprintf(stderr, "Integer out of range\n"); + return FALSE; + } + iptr = (int*)num; + *iptr = (int)tmp; + } + return TRUE; +} + +// the same as myatoll but for double +// There's no NAN & INF checking here (what if they would be needed?) +bool myatod(void *num, const char *str, argtype t){ + double tmp, *dptr; + float *fptr; + char *endptr; + assert(str); + tmp = strtod(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0') + return FALSE; + switch(t){ + case arg_double: + dptr = (double *) num; + *dptr = tmp; + break; + case arg_float: + default: + fptr = (float *) num; + *fptr = (float)tmp; + break; + } + return TRUE; +} + +/** + * Get index of current option in array options + * @param opt (i) - returning val of getopt_long + * @param options (i) - array of options + * @return index in array + */ +int get_optind(int opt, myoption *options){ + int oind; + myoption *opts = options; + assert(opts); + for(oind = 0; opts->name && opts->val != opt; oind++, opts++); + if(!opts->name || opts->val != opt) // no such parameter + showhelp(-1, options); + return oind; +} + +/** + * Parce command line arguments + * ! If arg is string, then value will be strdup'ed! + * + * @param argc (io) - address of argc of main(), return value of argc stay after `getopt` + * @param argv (io) - address of argv of main(), return pointer to argv stay after `getopt` + * BE CAREFUL! if you wanna use full argc & argv, save their original values before + * calling this function + * @param options (i) - array of `myoption` for arguments parcing + * + * @exit: in case of error this function show help & make `exit(-1)` + */ +void parceargs(int *argc, char ***argv, myoption *options){ + char *short_options, *soptr; + struct option *long_options, *loptr; + size_t optsize, i; + myoption *opts = options; + // check whether there is at least one options + assert(opts); + assert(opts[0].name); + // first we count how much values are in opts + for(optsize = 0; opts->name; optsize++, opts++); + // now we can allocate memory + short_options = calloc(optsize * 3 + 1, 1); // multiply by three for '::' in case of args in opts + long_options = calloc(optsize + 1, sizeof(struct option)); + opts = options; loptr = long_options; soptr = short_options; + // fill short/long parameters and make a simple checking + for(i = 0; i < optsize; i++, loptr++, opts++){ + // check + assert(opts->name); // check name + if(opts->has_arg){ + assert(opts->type != arg_none); // check error with arg type + assert(opts->argptr); // check pointer + } + if(opts->type != arg_none) // if there is a flag without arg, check its pointer + assert(opts->argptr); + // fill long_options + // don't do memcmp: what if there would be different alignment? + loptr->name = opts->name; + loptr->has_arg = opts->has_arg; + loptr->flag = opts->flag; + loptr->val = opts->val; + // fill short options if they are: + if(!opts->flag){ + *soptr++ = opts->val; + if(opts->has_arg) // add ':' if option has required argument + *soptr++ = ':'; + if(opts->has_arg == 2) // add '::' if option has optional argument + *soptr++ = ':'; + } + } + // now we have both long_options & short_options and can parse `getopt_long` + while(1){ + int opt; + int oindex = 0, optind = 0; // oindex - number of option in argv, optind - number in options[] + if((opt = getopt_long(*argc, *argv, short_options, long_options, &oindex)) == -1) break; + if(opt == '?'){ + opt = optopt; + optind = get_optind(opt, options); + if(options[optind].has_arg == 1) showhelp(optind, options); // need argument + } + else{ + if(opt == 0 || oindex > 0) optind = oindex; + else optind = get_optind(opt, options); + } + opts = &options[optind]; +DBG ("\n*******\noption %s (oindex = %d / optind = %d)", options[optind].name, oindex, optind); +if(optarg) DBG (" with arg %s", optarg); +DBG ("\n"); + if(opt == 0 && opts->has_arg == 0) continue; // only long option changing integer flag +DBG("opt = %c, arg type: ", opt); + // now check option + if(opts->has_arg == 1) assert(optarg); + bool result = TRUE; + // even if there is no argument, but argptr != NULL, think that optarg = "1" + if(!optarg) optarg = "1"; + switch(opts->type){ + default: + case arg_none: +DBG("none\n"); + break; + case arg_int: +DBG("integer\n"); + result = myatoll(opts->argptr, optarg, arg_int); + break; + case arg_longlong: +DBG("long long\n"); + result = myatoll(opts->argptr, optarg, arg_longlong); + break; + case arg_double: +DBG("double\n"); + result = myatod(opts->argptr, optarg, arg_double); + break; + case arg_float: +DBG("double\n"); + result = myatod(opts->argptr, optarg, arg_float); + break; + case arg_string: +DBG("string\n"); + result = (*((char **)opts->argptr) = strdup(optarg)); + break; + case arg_function: +DBG("function\n"); + result = ((argfn)opts->argptr)(optarg, optind); + break; + } + if(!result){ +DBG("OOOPS! Error in result\n"); + showhelp(optind, options); + } + } + *argc -= optind; + *argv += optind; +} + +/** + * Show help information based on myoption->help values + * @param oindex (i) - if non-negative, show only help by myoption[oindex].help + * @param options (i) - array of `myoption` + * + * @exit: run `exit(-1)` !!! + */ +void showhelp(int oindex, myoption *options){ + // ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext, + // but you can redefine it before `#include "parceargs.h"` + int max_opt_len = 0; // max len of options substring - for right indentation + const int bufsz = 255; + char buf[bufsz+1]; + myoption *opts = options; + assert(opts); + assert(opts[0].name); // check whether there is at least one options + if(oindex > -1){ // print only one message + opts = &options[oindex]; + printf(" "); + if(!opts->flag && isalpha(opts->val)) printf("-%c, ", opts->val); + printf("--%s", opts->name); + if(opts->has_arg == 1) printf("=arg"); + else if(opts->has_arg == 2) printf("[=arg]"); + printf(" %s\n", PRNT(opts->help)); + exit(-1); + } + // header, by default is just "progname\n" + printf("\n"); + if(strstr(helpstring, "%s")) // print progname + printf(helpstring, __progname); + else // only text + printf("%s", helpstring); + printf("\n"); + // count max_opt_len + do{ + int L = strlen(opts->name); + if(max_opt_len < L) max_opt_len = L; + }while((++opts)->name); + max_opt_len += 14; // format: '-S , --long[=arg]' - get addition 13 symbols + opts = options; + // Now print all help + do{ + int p = sprintf(buf, " "); // a little indent + if(!opts->flag && isalpha(opts->val)) // .val is short argument + p += snprintf(buf+p, bufsz-p, "-%c, ", opts->val); + p += snprintf(buf+p, bufsz-p, "--%s", opts->name); + if(opts->has_arg == 1) // required argument + p += snprintf(buf+p, bufsz-p, "=arg"); + else if(opts->has_arg == 2) // optional argument + p += snprintf(buf+p, bufsz-p, "[=arg]"); + assert(p < max_opt_len); // there would be magic if p >= max_opt_len + printf("%-*s%s\n", max_opt_len+1, buf, PRNT(opts->help)); // write options & at least 2 spaces after + }while((++opts)->name); + printf("\n\n"); + exit(-1); +} diff --git a/kicad-copy/parceargs.h b/kicad-copy/parceargs.h new file mode 100644 index 0000000..554803b --- /dev/null +++ b/kicad-copy/parceargs.h @@ -0,0 +1,105 @@ +/* + * parceargs.h - headers for parcing command line arguments + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __PARCEARGS_H__ +#define __PARCEARGS_H__ + +#include // bool +#include + +#ifndef TRUE + #define TRUE true +#endif + +#ifndef FALSE + #define FALSE false +#endif + +// macro for argptr +#define APTR(x) ((void*)x) + +// if argptr is a function: +typedef bool(*argfn)(void *arg, int N); + +/* + * type of getopt's argument + * WARNING! + * My function change value of flags by pointer, so if you want to use another type + * make a latter conversion, example: + * char charg; + * int iarg; + * myoption opts[] = { + * {"value", 1, NULL, 'v', arg_int, &iarg, "char val"}, ..., end_option}; + * ..(parce args).. + * charg = (char) iarg; + */ +typedef enum { + arg_none = 0, // no arg + arg_int, // integer + arg_longlong, // long long + arg_double, // double + arg_float, // float + arg_string, // char * + arg_function // parce_args will run function `bool (*fn)(char *optarg, int N)` +} argtype; + +/* + * Structure for getopt_long & help + * BE CAREFUL: .argptr is pointer to data or pointer to function, + * conversion depends on .type + * + * ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext, + * but you can redefine it before `#include "parceargs.h"` + * + * if arg is string, then value wil be strdup'ed like that: + * char *str; + * myoption opts[] = {{"string", 1, NULL, 's', arg_string, &str, "string val"}, ..., end_option}; + * *(opts[1].str) = strdup(optarg); + * in other cases argptr should be address of some variable (or pointer to allocated memory) + * + * NON-NULL argptr should be written inside macro APTR(argptr) or directly: (void*)argptr + * + * !!!LAST VALUE OF ARRAY SHOULD BE `end_option` or ZEROS !!! + * + */ +typedef struct{ + // these are from struct option: + const char *name; // long option's name + int has_arg; // 0 - no args, 1 - nesessary arg, 2 - optionally arg + int *flag; // NULL to return val, pointer to int - to set its value of val (function returns 0) + int val; // short opt name (if flag == NULL) or flag's value + // and these are mine: + argtype type; // type of argument + void *argptr; // pointer to variable to assign optarg value or function `bool (*fn)(char *optarg, int N)` + char *help; // help string which would be shown in function `showhelp` or NULL +} myoption; + +// last string of array (all zeros) +#define end_option {0,0,0,0,0,0,0} + + +extern const char *__progname; + +void showhelp(int oindex, myoption *options); +void parceargs(int *argc, char ***argv, myoption *options); +void change_helpstring(char *s); + +#endif // __PARCEARGS_H__ diff --git a/kicad-copy/usefull_macros.c b/kicad-copy/usefull_macros.c new file mode 100644 index 0000000..524c771 --- /dev/null +++ b/kicad-copy/usefull_macros.c @@ -0,0 +1,300 @@ +/* + * usefull_macros.h - a set of usefull functions: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include "usefull_macros.h" + +/******************************************************************************\ + * Coloured terminal +\******************************************************************************/ +int globErr = 0; // errno for WARN/ERR + +// pointers to coloured output printf +int (*red)(const char *fmt, ...); +int (*green)(const char *fmt, ...); +int (*_WARN)(const char *fmt, ...); + +/* + * format red / green messages + * name: r_pr_, g_pr_ + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(RED); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +int g_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(GREEN); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +/* + * print red error/warning messages (if output is a tty) + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_WARN(const char *fmt, ...){ + va_list ar; int i = 1; + fprintf(stderr, RED); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = vfprintf(stderr, fmt, ar); + va_end(ar); + i++; + fprintf(stderr, OLDCOLOR "\n"); + return i; +} + +const char stars[] = "****************************************"; +/* + * notty variants of coloured printf + * name: s_WARN, r_pr_notty + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int s_WARN(const char *fmt, ...){ + va_list ar; int i; + i = fprintf(stderr, "\n%s\n", stars); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = +vfprintf(stderr, fmt, ar); + va_end(ar); + i += fprintf(stderr, "\n%s\n", stars); + i += fprintf(stderr, "\n"); + return i; +} +int r_pr_notty(const char *fmt, ...){ + va_list ar; int i; + i = printf("\n%s\n", stars); + va_start(ar, fmt); + i += vprintf(fmt, ar); + va_end(ar); + i += printf("\n%s\n", stars); + return i; +} + +/** + * Run this function in the beginning of main() to setup locale & coloured output + */ +void initial_setup(){ + // setup coloured output + if(isatty(STDOUT_FILENO)){ // make color output in tty + red = r_pr_; green = g_pr_; + }else{ // no colors in case of pipe + red = r_pr_notty; green = printf; + } + if(isatty(STDERR_FILENO)) _WARN = r_WARN; + else _WARN = s_WARN; + // Setup locale + setlocale(LC_ALL, ""); + setlocale(LC_NUMERIC, "C"); + bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR); + textdomain(GETTEXT_PACKAGE); +} + +/******************************************************************************\ + * Memory +\******************************************************************************/ +/* + * safe memory allocation for macro ALLOC + * @param N - number of elements to allocate + * @param S - size of single element (typically sizeof) + * @return pointer to allocated memory area + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) ERR("malloc"); + //assert(p); + return p; +} +/* +#include +#include +#include +#include +*/ +/** + * Mmap file to a memory area + * + * @param filename (i) - name of file to mmap + * @return stuct with mmap'ed file or die + */ +mmapbuf *My_mmap(char *filename){ + int fd; + char *ptr; + size_t Mlen; + struct stat statbuf; + if(!filename) ERRX(_("No filename given!")); + if((fd = open(filename, O_RDONLY)) < 0) + ERR(_("Can't open %s for reading"), filename); + if(fstat (fd, &statbuf) < 0) + ERR(_("Can't stat %s"), filename); + Mlen = statbuf.st_size; + if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + ERR(_("Mmap error for input")); + if(close(fd)) ERR(_("Can't close mmap'ed file")); + mmapbuf *ret = MALLOC(mmapbuf, 1); + ret->data = ptr; + ret->len = Mlen; + return ret; +} + +void My_munmap(mmapbuf *b){ + if(munmap(b->data, b->len)) + ERR(_("Can't munmap")); + FREE(b); +} + + +/******************************************************************************\ + * Terminal in no-echo mode +\******************************************************************************/ +struct termios oldt, newt; // terminal flags +// run on exit: +/* +void quit(int sig){ + //... + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state + //... +} +*/ +// initial setup: +void setup_con(){ + tcgetattr(STDIN_FILENO, &oldt); + newt = oldt; + newt.c_lflag &= ~(ICANON | ECHO); + if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0){ + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); + exit(-2); //quit? + } +} + +/** + * Read character from console without echo + * @return char readed + */ +int read_console(){ + int rb; + struct timeval tv; + int retval; + fd_set rfds; + FD_ZERO(&rfds); + FD_SET(STDIN_FILENO, &rfds); + tv.tv_sec = 0; tv.tv_usec = 10000; + retval = select(1, &rfds, NULL, NULL, &tv); + if(!retval) rb = 0; + else { + if(FD_ISSET(STDIN_FILENO, &rfds)) rb = getchar(); + else rb = 0; + } + return rb; +} + +/** + * getchar() without echo + * wait until at least one character pressed + * @return character readed + */ +int mygetchar(){ // аналог getchar() без необходимости жать Enter + int ret; + do ret = read_console(); + while(ret == 0); + return ret; +} + + +/******************************************************************************\ + * TTY with select() +\******************************************************************************/ +struct termio oldtty, tty; // TTY flags +char *comdev; // TTY device name +int comfd; // TTY fd +// run on exit: +/* +void quit(int ex_stat){ + ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state + close(comfd); + //... +} +*/ +#ifndef BAUD_RATE +#define BAUD_RATE B9600 +#endif +// init: +void tty_init(){ + printf("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + fprintf(stderr,"Can't use port %s\n",comdev); + ioctl(comfd, TCSANOW, &oldtty); // return TTY to previous state + close(comfd); + exit(1); // quit? + } + printf(" OK\nGet current settings...\n"); + if(ioctl(comfd,TCGETA,&oldtty) < 0) exit(-1); // Get settings + tty = oldtty; + tty.c_lflag = 0; // ~(ICANON | ECHO | ECHOE | ISIG) + tty.c_oflag = 0; + tty.c_cflag = BAUD_RATE|CS8|CREAD|CLOCAL; // 9.6k, 8N1, RW, ignore line ctrl + tty.c_cc[VMIN] = 0; // non-canonical mode + tty.c_cc[VTIME] = 5; + if(ioctl(comfd,TCSETA,&tty) < 0) exit(-1); // set new mode + printf(" OK\n"); +} +/** + * Read data from TTY + * @param buff (o) - buffer for data read + * @param length - buffer len + * @return amount of readed bytes + */ +size_t read_tty(uint8_t *buff, size_t length){ + ssize_t L = 0; + fd_set rfds; + struct timeval tv; + int retval; + FD_ZERO(&rfds); + FD_SET(comfd, &rfds); + tv.tv_sec = 0; tv.tv_usec = 50000; // wait for 50ms + retval = select(comfd + 1, &rfds, NULL, NULL, &tv); + if (!retval) return 0; + if(FD_ISSET(comfd, &rfds)){ + if((L = read(comfd, buff, length)) < 1) return 0; + } + return (size_t)L; +} diff --git a/kicad-copy/usefull_macros.h b/kicad-copy/usefull_macros.h new file mode 100644 index 0000000..43a6192 --- /dev/null +++ b/kicad-copy/usefull_macros.h @@ -0,0 +1,106 @@ +/* + * usefull_macros.h - a set of usefull macros: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __USEFULL_MACROS_H__ +#define __USEFULL_MACROS_H__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* + * GETTEXT + */ +#define _(String) gettext(String) +#define gettext_noop(String) String +#define N_(String) gettext_noop(String) + +// unused arguments with -Wall -Werror +#define _U_ __attribute__((__unused__)) + +/* + * Coloured messages output + */ +#define RED "\033[1;31;40m" +#define GREEN "\033[1;32;40m" +#define OLDCOLOR "\033[0;0;0m" + +/* + * ERROR/WARNING messages + */ +extern int globErr; +#define ERR(...) do{globErr=errno; _WARN(__VA_ARGS__); exit(-1);}while(0) +#define ERRX(...) do{globErr=0; _WARN(__VA_ARGS__); exit(-1);}while(0) +#define WARN(...) do{globErr=errno; _WARN(__VA_ARGS__);}while(0) +#define WARNX(...) do{globErr=0; _WARN(__VA_ARGS__);}while(0) + +/* + * print function name, debug messages + * debug mode, -DEBUG + */ +#ifdef EBUG + #define FNAME() fprintf(stderr, "\n%s (%s, line %d)\n", __func__, __FILE__, __LINE__) + #define DBG(...) do{fprintf(stderr, "%s (%s, line %d): ", __func__, __FILE__, __LINE__); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n");} while(0) +#else + #define FNAME() do{}while(0) + #define DBG(...) do{}while(0) +#endif //EBUG + +/* + * Memory allocation + */ +#define ALLOC(type, var, size) type * var = ((type *)my_alloc(size, sizeof(type))) +#define MALLOC(type, size) ((type *)my_alloc(size, sizeof(type))) +#define FREE(ptr) do{free(ptr); ptr = NULL;}while(0) + +// functions for color output in tty & no-color in pipes +extern int (*red)(const char *fmt, ...); +extern int (*_WARN)(const char *fmt, ...); +extern int (*green)(const char *fmt, ...); +void * my_alloc(size_t N, size_t S); +void initial_setup(); + +// mmap file +typedef struct{ + char *data; + size_t len; +} mmapbuf; +mmapbuf *My_mmap(char *filename); +void My_munmap(mmapbuf *b); + +#endif // __USEFULL_MACROS_H__ diff --git a/simple_gettext/1.c b/simple_gettext/1.c new file mode 100644 index 0000000..f19a854 --- /dev/null +++ b/simple_gettext/1.c @@ -0,0 +1,32 @@ +#include +#include +#include + +#define _(String) gettext(String) +#define gettext_noop(String) String +#define N_(String) gettext_noop(String) + +int main(int argc, char **argv){ + setlocale(LC_ALL, ""); + setlocale(LC_NUMERIC, "C"); + bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR); + /* + * To run in GUI bullshit (like qt or gtk) you must do: + * bind_textdomain_codeset(PACKAGE, "UTF8"); + */ + textdomain(GETTEXT_PACKAGE); + /* + * In case when you need both GUI and CLI, you should + * do this: + * char*old = bind_textdomain_codeset (PACKAGE, ""); + * before printf calling + * To restore GUI hrunicode do + * bind_textdomain_codeset (PACKAGE, old); + * at the end of printf's + */ + /// \n + printf(_("Translated by gettext\n")); + /// \n + printf(_("This works fine in any charset\n")); + return 0; +} diff --git a/simple_gettext/Makefile b/simple_gettext/Makefile new file mode 100644 index 0000000..c89277c --- /dev/null +++ b/simple_gettext/Makefile @@ -0,0 +1,31 @@ +PROGRAM = test +LDFLAGS = +DEFINES = -D_XOPEN_SOURCE=501 -DGETTEXT_PACKAGE=\"$(PROGRAM)\" -DLOCALEDIR=\".\" +SRCS = 1.c +CC = gcc +CXX = gcc +CFLAGS = -Wall -Werror $(DEFINES) +OBJS = $(SRCS:.c=.o) +PO_FILE = $(PROGRAM).po +RU_FILE = $(PROGRAM)_ru.po +LDIR = "ru/LC_MESSAGES" +MO_FILE = $(LDIR)/$(PROGRAM).mo + +all : $(PROGRAM) $(MO_FILE) + +$(PROGRAM) : $(OBJS) + $(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM) + +$(PO_FILE): $(SRCS) + xgettext -D. --from-code=koi8-r $(SRCS) -c -k_ -kN_ -o $(PO_FILE) + sed -i 's/charset=.*\\n/charset=koi8-r\\n/' $(PO_FILE) && enconv $(PO_FILE) +$(RU_FILE): $(PO_FILE) + [ -e $(RU_FILE) ] && msgmerge -Uis $(RU_FILE) $(PO_FILE) || cp $(PO_FILE) $(RU_FILE) +$(MO_FILE): $(LDIR) $(RU_FILE) + msgfmt $(RU_FILE) -o $(MO_FILE) +$(LDIR): + [ -e $(LDIR) ] || mkdir -p $(LDIR) +clean: + /bin/rm -f *.o *~ +depend: + $(CXX) -MM $(CXX.SRCS) diff --git a/simple_gettext/test_ru.po b/simple_gettext/test_ru.po new file mode 100644 index 0000000..7d6b857 --- /dev/null +++ b/simple_gettext/test_ru.po @@ -0,0 +1,29 @@ +# SOME DESCRIPTIVE TITLE. +# Copyright (C) YEAR THE PACKAGE'S COPYRIGHT HOLDER +# This file is distributed under the same license as the PACKAGE package. +# FIRST AUTHOR , YEAR. +# +#, fuzzy +msgid "" +msgstr "Project-Id-Version: PACKAGE VERSION\n" + "Report-Msgid-Bugs-To: \n" + "POT-Creation-Date: 2014-04-24 18:13+0400\n" + "PO-Revision-Date: YEAR-MO-DA HO:MI+ZONE\n" + "Last-Translator: FULL NAME \n" + "Language-Team: LANGUAGE \n" + "Language: \n" + "MIME-Version: 1.0\n" + "Content-Type: text/plain; charset=koi8-r\n" + "Content-Transfer-Encoding: 8bit\n" + +#. / \n +#: 1.c:17 +#, c-format +msgid "This works fine in any charset\n" +msgstr " \n" + +#. / \n +#: 1.c:15 +#, c-format +msgid "Translated by gettext\n" +msgstr " \n" diff --git a/simple_list.c b/simple_list.c new file mode 100644 index 0000000..6e237a1 --- /dev/null +++ b/simple_list.c @@ -0,0 +1,112 @@ +/* + * simple_list.c - simple one-direction list + * + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include +#include +#include +#include + +#include "simple_list.h" + +#define MALLOC alloc_simple_list +/** + * Memory allocation with control + * @param N - number of elements + * @param S - size of one element + * @return allocated memory or die + */ +static void *alloc_simple_list(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) err(1, "malloc"); + return p; +} + +void (*listdata_free)(listdata node_data) = NULL; // external function to free(listdata) + +#define FREE(arg) do{free(arg); arg = NULL;}while(0) + +/** + * add element v to list with root root (also this can be last element) + * @param root (io) - pointer to root (or last element) of list or NULL + * if *root == NULL, just created node will be placed there + * @param v - data inserted + * @return pointer to created node + */ +List *list_add(List **root, listdata v){ + List *node, *last; + if((node = (List*) MALLOC(1, sizeof(List))) == 0) return NULL; // allocation error + node->data = v; // insert data + if(root){ + if(*root){ // there was root node - search last + last = *root; + while(last->next) last = last->next; + last->next = node; // insert pointer to new node into last element in list + } + if(!*root) *root = node; + } + return node; +} + +/** + * set listdata_free() + * @param fn - function that freec memory of (node) + */ +void listfree_function(void (*fn)(listdata node)){ + listdata_free = fn; +} + +/** + * remove all nodes in list + * @param root - pointer to root node + */ +void list_free(List **root){ + List *node = *root, *next; + if(!root || !*root) return; + do{ + next = node->next; + if(listdata_free) + listdata_free(node->data); + free(node); + node = next; + }while(node); + *root = NULL; +} + +#ifdef STANDALONE +int main(void) { + List *tp = NULL, *root_p = NULL; + int i, ins[] = {4,2,6,1,3,4,7}; + for(i = 0; i < 7; i++){ + if(!(tp = list_add(&tp, ins[i]))) + err(1, "Malloc error"); // can't insert + if(!root_p) root_p = tp; + } + tp = root_p; + i = 0; + do{ + printf("element %d = %d\n", i++, tp->data); + tp = tp->next; + }while(tp); + list_free(&root_p); + return 0; +} +#endif diff --git a/simple_list.h b/simple_list.h new file mode 100644 index 0000000..e8cdc68 --- /dev/null +++ b/simple_list.h @@ -0,0 +1,44 @@ +/* + * simple_list.h - header file for simple list support + * TO USE IT you must define the type of data as + * typedef your_type listdata + * or at compiling time + * -Dlistdata=your_type + * or by changing this file + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ +#pragma once +#ifndef __SIMPLE_LIST_H__ +#define __SIMPLE_LIST_H__ + +#ifdef STANDALONE +typedef int listdata; +#endif + +typedef struct list_{ + listdata data; + struct list_ *next; +} List; + +// add element v to list with root root (also this can be last element) +List *list_add(List **root, listdata v); +// set listdata_free() +void listfree_function(void (*fn)(listdata node)); + +#endif // __SIMPLE_LIST_H__ diff --git a/usefull_macros.c b/usefull_macros.c new file mode 100644 index 0000000..524c771 --- /dev/null +++ b/usefull_macros.c @@ -0,0 +1,300 @@ +/* + * usefull_macros.h - a set of usefull functions: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#include "usefull_macros.h" + +/******************************************************************************\ + * Coloured terminal +\******************************************************************************/ +int globErr = 0; // errno for WARN/ERR + +// pointers to coloured output printf +int (*red)(const char *fmt, ...); +int (*green)(const char *fmt, ...); +int (*_WARN)(const char *fmt, ...); + +/* + * format red / green messages + * name: r_pr_, g_pr_ + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(RED); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +int g_pr_(const char *fmt, ...){ + va_list ar; int i; + printf(GREEN); + va_start(ar, fmt); + i = vprintf(fmt, ar); + va_end(ar); + printf(OLDCOLOR); + return i; +} +/* + * print red error/warning messages (if output is a tty) + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int r_WARN(const char *fmt, ...){ + va_list ar; int i = 1; + fprintf(stderr, RED); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = vfprintf(stderr, fmt, ar); + va_end(ar); + i++; + fprintf(stderr, OLDCOLOR "\n"); + return i; +} + +const char stars[] = "****************************************"; +/* + * notty variants of coloured printf + * name: s_WARN, r_pr_notty + * @param fmt ... - printf-like format + * @return number of printed symbols + */ +int s_WARN(const char *fmt, ...){ + va_list ar; int i; + i = fprintf(stderr, "\n%s\n", stars); + va_start(ar, fmt); + if(globErr){ + errno = globErr; + vwarn(fmt, ar); + errno = 0; + globErr = 0; + }else + i = +vfprintf(stderr, fmt, ar); + va_end(ar); + i += fprintf(stderr, "\n%s\n", stars); + i += fprintf(stderr, "\n"); + return i; +} +int r_pr_notty(const char *fmt, ...){ + va_list ar; int i; + i = printf("\n%s\n", stars); + va_start(ar, fmt); + i += vprintf(fmt, ar); + va_end(ar); + i += printf("\n%s\n", stars); + return i; +} + +/** + * Run this function in the beginning of main() to setup locale & coloured output + */ +void initial_setup(){ + // setup coloured output + if(isatty(STDOUT_FILENO)){ // make color output in tty + red = r_pr_; green = g_pr_; + }else{ // no colors in case of pipe + red = r_pr_notty; green = printf; + } + if(isatty(STDERR_FILENO)) _WARN = r_WARN; + else _WARN = s_WARN; + // Setup locale + setlocale(LC_ALL, ""); + setlocale(LC_NUMERIC, "C"); + bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR); + textdomain(GETTEXT_PACKAGE); +} + +/******************************************************************************\ + * Memory +\******************************************************************************/ +/* + * safe memory allocation for macro ALLOC + * @param N - number of elements to allocate + * @param S - size of single element (typically sizeof) + * @return pointer to allocated memory area + */ +void *my_alloc(size_t N, size_t S){ + void *p = calloc(N, S); + if(!p) ERR("malloc"); + //assert(p); + return p; +} +/* +#include +#include +#include +#include +*/ +/** + * Mmap file to a memory area + * + * @param filename (i) - name of file to mmap + * @return stuct with mmap'ed file or die + */ +mmapbuf *My_mmap(char *filename){ + int fd; + char *ptr; + size_t Mlen; + struct stat statbuf; + if(!filename) ERRX(_("No filename given!")); + if((fd = open(filename, O_RDONLY)) < 0) + ERR(_("Can't open %s for reading"), filename); + if(fstat (fd, &statbuf) < 0) + ERR(_("Can't stat %s"), filename); + Mlen = statbuf.st_size; + if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + ERR(_("Mmap error for input")); + if(close(fd)) ERR(_("Can't close mmap'ed file")); + mmapbuf *ret = MALLOC(mmapbuf, 1); + ret->data = ptr; + ret->len = Mlen; + return ret; +} + +void My_munmap(mmapbuf *b){ + if(munmap(b->data, b->len)) + ERR(_("Can't munmap")); + FREE(b); +} + + +/******************************************************************************\ + * Terminal in no-echo mode +\******************************************************************************/ +struct termios oldt, newt; // terminal flags +// run on exit: +/* +void quit(int sig){ + //... + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state + //... +} +*/ +// initial setup: +void setup_con(){ + tcgetattr(STDIN_FILENO, &oldt); + newt = oldt; + newt.c_lflag &= ~(ICANON | ECHO); + if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0){ + tcsetattr(STDIN_FILENO, TCSANOW, &oldt); + exit(-2); //quit? + } +} + +/** + * Read character from console without echo + * @return char readed + */ +int read_console(){ + int rb; + struct timeval tv; + int retval; + fd_set rfds; + FD_ZERO(&rfds); + FD_SET(STDIN_FILENO, &rfds); + tv.tv_sec = 0; tv.tv_usec = 10000; + retval = select(1, &rfds, NULL, NULL, &tv); + if(!retval) rb = 0; + else { + if(FD_ISSET(STDIN_FILENO, &rfds)) rb = getchar(); + else rb = 0; + } + return rb; +} + +/** + * getchar() without echo + * wait until at least one character pressed + * @return character readed + */ +int mygetchar(){ // аналог getchar() без необходимости жать Enter + int ret; + do ret = read_console(); + while(ret == 0); + return ret; +} + + +/******************************************************************************\ + * TTY with select() +\******************************************************************************/ +struct termio oldtty, tty; // TTY flags +char *comdev; // TTY device name +int comfd; // TTY fd +// run on exit: +/* +void quit(int ex_stat){ + ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state + close(comfd); + //... +} +*/ +#ifndef BAUD_RATE +#define BAUD_RATE B9600 +#endif +// init: +void tty_init(){ + printf("\nOpen port...\n"); + if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){ + fprintf(stderr,"Can't use port %s\n",comdev); + ioctl(comfd, TCSANOW, &oldtty); // return TTY to previous state + close(comfd); + exit(1); // quit? + } + printf(" OK\nGet current settings...\n"); + if(ioctl(comfd,TCGETA,&oldtty) < 0) exit(-1); // Get settings + tty = oldtty; + tty.c_lflag = 0; // ~(ICANON | ECHO | ECHOE | ISIG) + tty.c_oflag = 0; + tty.c_cflag = BAUD_RATE|CS8|CREAD|CLOCAL; // 9.6k, 8N1, RW, ignore line ctrl + tty.c_cc[VMIN] = 0; // non-canonical mode + tty.c_cc[VTIME] = 5; + if(ioctl(comfd,TCSETA,&tty) < 0) exit(-1); // set new mode + printf(" OK\n"); +} +/** + * Read data from TTY + * @param buff (o) - buffer for data read + * @param length - buffer len + * @return amount of readed bytes + */ +size_t read_tty(uint8_t *buff, size_t length){ + ssize_t L = 0; + fd_set rfds; + struct timeval tv; + int retval; + FD_ZERO(&rfds); + FD_SET(comfd, &rfds); + tv.tv_sec = 0; tv.tv_usec = 50000; // wait for 50ms + retval = select(comfd + 1, &rfds, NULL, NULL, &tv); + if (!retval) return 0; + if(FD_ISSET(comfd, &rfds)){ + if((L = read(comfd, buff, length)) < 1) return 0; + } + return (size_t)L; +} diff --git a/usefull_macros.h b/usefull_macros.h new file mode 100644 index 0000000..43a6192 --- /dev/null +++ b/usefull_macros.h @@ -0,0 +1,106 @@ +/* + * usefull_macros.h - a set of usefull macros: memory, color etc + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +#pragma once +#ifndef __USEFULL_MACROS_H__ +#define __USEFULL_MACROS_H__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* + * GETTEXT + */ +#define _(String) gettext(String) +#define gettext_noop(String) String +#define N_(String) gettext_noop(String) + +// unused arguments with -Wall -Werror +#define _U_ __attribute__((__unused__)) + +/* + * Coloured messages output + */ +#define RED "\033[1;31;40m" +#define GREEN "\033[1;32;40m" +#define OLDCOLOR "\033[0;0;0m" + +/* + * ERROR/WARNING messages + */ +extern int globErr; +#define ERR(...) do{globErr=errno; _WARN(__VA_ARGS__); exit(-1);}while(0) +#define ERRX(...) do{globErr=0; _WARN(__VA_ARGS__); exit(-1);}while(0) +#define WARN(...) do{globErr=errno; _WARN(__VA_ARGS__);}while(0) +#define WARNX(...) do{globErr=0; _WARN(__VA_ARGS__);}while(0) + +/* + * print function name, debug messages + * debug mode, -DEBUG + */ +#ifdef EBUG + #define FNAME() fprintf(stderr, "\n%s (%s, line %d)\n", __func__, __FILE__, __LINE__) + #define DBG(...) do{fprintf(stderr, "%s (%s, line %d): ", __func__, __FILE__, __LINE__); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n");} while(0) +#else + #define FNAME() do{}while(0) + #define DBG(...) do{}while(0) +#endif //EBUG + +/* + * Memory allocation + */ +#define ALLOC(type, var, size) type * var = ((type *)my_alloc(size, sizeof(type))) +#define MALLOC(type, size) ((type *)my_alloc(size, sizeof(type))) +#define FREE(ptr) do{free(ptr); ptr = NULL;}while(0) + +// functions for color output in tty & no-color in pipes +extern int (*red)(const char *fmt, ...); +extern int (*_WARN)(const char *fmt, ...); +extern int (*green)(const char *fmt, ...); +void * my_alloc(size_t N, size_t S); +void initial_setup(); + +// mmap file +typedef struct{ + char *data; + size_t len; +} mmapbuf; +mmapbuf *My_mmap(char *filename); +void My_munmap(mmapbuf *b); + +#endif // __USEFULL_MACROS_H__