From 8602fda18766eb009b0bf66cb777e8a6a3c96886 Mon Sep 17 00:00:00 2001 From: Eddy Date: Wed, 19 Mar 2014 00:57:05 +0400 Subject: [PATCH] copy --- CMakeLists.txt | 34 + COPYING | 2 + ChangeLog | 0 INSTALL | 4 + NEWS | 1 + RUN | 4 + mask/README.fmt | 17 + mask/holes.json | 264 ++++++++ mask/json.c | 205 ++++++ mask/make_mask.m | 34 + src/CMakeLists.txt | 107 ++++ src/CPU.c | 113 ++++ src/CUDA.cu | 807 ++++++++++++++++++++++++ src/README | 13 + src/cmdlnopts.c | 299 +++++++++ src/diaphragm.c | 298 +++++++++ src/include/cmdlnopts.h | 48 ++ src/include/cutil_math.h | 772 +++++++++++++++++++++++ src/include/diaphragm.h | 31 + src/include/mkHartmann.h | 187 ++++++ src/include/parceargs.h | 104 +++ src/include/saveimg.h | 42 ++ src/include/wrapper.h | 144 +++++ src/locale/ru/LC_MESSAGES/mkHartmann.mo | Bin 0 -> 1118 bytes src/locale/ru/messages.po | 304 +++++++++ src/locale/ru/ru.po | 305 +++++++++ src/locale/ru/ru.po~ | 301 +++++++++ src/locale/ru/update | 2 + src/mkHartmann.c | 446 +++++++++++++ src/parceargs.c | 316 ++++++++++ src/saveimg.c | 423 +++++++++++++ src/test/capability.cu | 25 + src/wrapper.c | 373 +++++++++++ 33 files changed, 6025 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 COPYING create mode 100644 ChangeLog create mode 100644 INSTALL create mode 100644 NEWS create mode 100755 RUN create mode 100644 mask/README.fmt create mode 100644 mask/holes.json create mode 100644 mask/json.c create mode 100644 mask/make_mask.m create mode 100644 src/CMakeLists.txt create mode 100644 src/CPU.c create mode 100644 src/CUDA.cu create mode 100644 src/README create mode 100644 src/cmdlnopts.c create mode 100644 src/diaphragm.c create mode 100644 src/include/cmdlnopts.h create mode 100644 src/include/cutil_math.h create mode 100644 src/include/diaphragm.h create mode 100644 src/include/mkHartmann.h create mode 100644 src/include/parceargs.h create mode 100644 src/include/saveimg.h create mode 100644 src/include/wrapper.h create mode 100644 src/locale/ru/LC_MESSAGES/mkHartmann.mo create mode 100644 src/locale/ru/messages.po create mode 100644 src/locale/ru/ru.po create mode 100644 src/locale/ru/ru.po~ create mode 100755 src/locale/ru/update create mode 100644 src/mkHartmann.c create mode 100644 src/parceargs.c create mode 100644 src/saveimg.c create mode 100644 src/test/capability.cu create mode 100644 src/wrapper.c diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..041c00a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,34 @@ +cmake_minimum_required(VERSION 2.8) +set(PROJ mkHartmann) +set(PROJ_VER "0.0.1") +project(${PROJ}) +set(CMAKE_COLOR_MAKEFILE ON) +#set(DEBUG 1) # Comment out this line in release mode +if(NOT DEFINED NO_CUDA) + message("Try to use CUDA") + find_package(CUDA) + if(CUDA_FOUND) + add_definitions(-DCUDA_FOUND) + endif() +endif() +if(NOT DEFINED PROCESSOR_COUNT) + set(PROCESSOR_COUNT 2) # by default 2 cores + set(cpuinfo_file "/proc/cpuinfo") + if(EXISTS "${cpuinfo_file}") + file(STRINGS "${cpuinfo_file}" procs REGEX "^processor.: [0-9]+$") + list(LENGTH procs PROCESSOR_COUNT) + endif() +endif() +add_definitions(-DTHREAD_NUMBER=${PROCESSOR_COUNT}) +message("In multithreaded operations will use ${PROCESSOR_COUNT} threads") + +# fix invalid install paths +if(DEFINED CMAKE_INSTALL_PREFIX AND CMAKE_INSTALL_PREFIX MATCHES "/usr/local") + set(CMAKE_INSTALL_PREFIX "/usr") +endif() +# change path by user's wish +if(NOT DEFINED LOCALEDIR) + set(LOCALEDIR "${CMAKE_INSTALL_PREFIX}/share/locale") +endif() +set(CMAKE_INSTALL_LOCALEDIR "${LOCALEDIR}/ru/LC_MESSAGES") +subdirs(src) diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..256a913 --- /dev/null +++ b/COPYING @@ -0,0 +1,2 @@ +eddy@sao.ru + diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 0000000..e69de29 diff --git a/INSTALL b/INSTALL new file mode 100644 index 0000000..b239d06 --- /dev/null +++ b/INSTALL @@ -0,0 +1,4 @@ +1. ./RUN +2. cd tmp +3. make +[4. make install] diff --git a/NEWS b/NEWS new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/NEWS @@ -0,0 +1 @@ + diff --git a/RUN b/RUN new file mode 100755 index 0000000..31ff689 --- /dev/null +++ b/RUN @@ -0,0 +1,4 @@ +#!/bin/bash +[ -d "tmp" ] && rm -rf tmp +mkdir tmp && cd tmp && cmake .. -DCMAKE_INSTALL_PREFIX=/home/eddy -DLOCALEDIR=/home/eddy/.local/share/locale -DDEBUG=1 +ln -s ../mask/holes.json diff --git a/mask/README.fmt b/mask/README.fmt new file mode 100644 index 0000000..b5c5a32 --- /dev/null +++ b/mask/README.fmt @@ -0,0 +1,17 @@ +JSON format: + +global object MUST consist field "Z" or "maskz" with Z-coordinate of mask + +it CAN contain global fields (common to all holes) such as in holes array +its field bbox (if exists) HAVE TO BE GLOBAL convex bbox of all holes! + +it MUST contain at least one hole object (or hole array) + +fields of "hole": + "shape" - shape of hole (could be "square" and "round" or "ellipse") + "radius" - scalar or array[2] with radius parameter + "center" - array[2] with coordinates of hole's center + "bbox" - array[4] with bounding box of a hole + + independently from shape any hole can consist both fields "radius" and "center" or single field "bbox" + \ No newline at end of file diff --git a/mask/holes.json b/mask/holes.json new file mode 100644 index 0000000..526f124 --- /dev/null +++ b/mask/holes.json @@ -0,0 +1,264 @@ +{ + "maskz": 20.0, + "shape": "round", "radius": 0.007500, + "holes": [ + { "ring": 0, "number": 0, "center": [ 0.1742, 0.0172 ] }, + { "ring": 0, "number": 1, "center": [ 0.1675, 0.0508 ] }, + { "ring": 0, "number": 2, "center": [ 0.1543, 0.0825 ] }, + { "ring": 0, "number": 3, "center": [ 0.1353, 0.1110 ] }, + { "ring": 0, "number": 4, "center": [ 0.1110, 0.1353 ] }, + { "ring": 0, "number": 5, "center": [ 0.0825, 0.1543 ] }, + { "ring": 0, "number": 6, "center": [ 0.0508, 0.1675 ] }, + { "ring": 0, "number": 7, "center": [ 0.0172, 0.1742 ] }, + { "ring": 0, "number": 8, "center": [ -0.0172, 0.1742 ] }, + { "ring": 0, "number": 9, "center": [ -0.0508, 0.1675 ] }, + { "ring": 0, "number": 10, "center": [ -0.0825, 0.1543 ] }, + { "ring": 0, "number": 11, "center": [ -0.1110, 0.1353 ] }, + { "ring": 0, "number": 12, "center": [ -0.1353, 0.1110 ] }, + { "ring": 0, "number": 13, "center": [ -0.1543, 0.0825 ] }, + { "ring": 0, "number": 14, "center": [ -0.1675, 0.0508 ] }, + { "ring": 0, "number": 15, "center": [ -0.1742, 0.0172 ] }, + { "ring": 0, "number": 16, "center": [ -0.1742, -0.0172 ] }, + { "ring": 0, "number": 17, "center": [ -0.1675, -0.0508 ] }, + { "ring": 0, "number": 18, "center": [ -0.1543, -0.0825 ] }, + { "ring": 0, "number": 19, "center": [ -0.1353, -0.1110 ] }, + { "ring": 0, "number": 20, "center": [ -0.1110, -0.1353 ] }, + { "ring": 0, "number": 21, "center": [ -0.0825, -0.1543 ] }, + { "ring": 0, "number": 22, "center": [ -0.0508, -0.1675 ] }, + { "ring": 0, "number": 23, "center": [ -0.0172, -0.1742 ] }, + { "ring": 0, "number": 24, "center": [ 0.0172, -0.1742 ] }, + { "ring": 0, "number": 25, "center": [ 0.0508, -0.1675 ] }, + { "ring": 0, "number": 26, "center": [ 0.0825, -0.1543 ] }, + { "ring": 0, "number": 27, "center": [ 0.1110, -0.1353 ] }, + { "ring": 0, "number": 28, "center": [ 0.1353, -0.1110 ] }, + { "ring": 0, "number": 29, "center": [ 0.1543, -0.0825 ] }, + { "ring": 0, "number": 30, "center": [ 0.1675, -0.0508 ] }, + { "ring": 0, "number": 31, "center": [ 0.1742, -0.0172 ] }, + { "ring": 1, "number": 0, "center": [ 0.2458, 0.0242 ] }, + { "ring": 1, "number": 1, "center": [ 0.2364, 0.0717 ] }, + { "ring": 1, "number": 2, "center": [ 0.2178, 0.1164 ] }, + { "ring": 1, "number": 3, "center": [ 0.1909, 0.1567 ] }, + { "ring": 1, "number": 4, "center": [ 0.1567, 0.1909 ] }, + { "ring": 1, "number": 5, "center": [ 0.1164, 0.2178 ] }, + { "ring": 1, "number": 6, "center": [ 0.0717, 0.2364 ] }, + { "ring": 1, "number": 7, "center": [ 0.0242, 0.2458 ] }, + { "ring": 1, "number": 8, "center": [ -0.0242, 0.2458 ] }, + { "ring": 1, "number": 9, "center": [ -0.0717, 0.2364 ] }, + { "ring": 1, "number": 10, "center": [ -0.1164, 0.2178 ] }, + { "ring": 1, "number": 11, "center": [ -0.1567, 0.1909 ] }, + { "ring": 1, "number": 12, "center": [ -0.1909, 0.1567 ] }, + { "ring": 1, "number": 13, "center": [ -0.2178, 0.1164 ] }, + { "ring": 1, "number": 14, "center": [ -0.2364, 0.0717 ] }, + { "ring": 1, "number": 15, "center": [ -0.2458, 0.0242 ] }, + { "ring": 1, "number": 16, "center": [ -0.2458, -0.0242 ] }, + { "ring": 1, "number": 17, "center": [ -0.2364, -0.0717 ] }, + { "ring": 1, "number": 18, "center": [ -0.2178, -0.1164 ] }, + { "ring": 1, "number": 19, "center": [ -0.1909, -0.1567 ] }, + { "ring": 1, "number": 20, "center": [ -0.1567, -0.1909 ] }, + { "ring": 1, "number": 21, "center": [ -0.1164, -0.2178 ] }, + { "ring": 1, "number": 22, "center": [ -0.0717, -0.2364 ] }, + { "ring": 1, "number": 23, "center": [ -0.0242, -0.2458 ] }, + { "ring": 1, "number": 24, "center": [ 0.0242, -0.2458 ] }, + { "ring": 1, "number": 25, "center": [ 0.0717, -0.2364 ] }, + { "ring": 1, "number": 26, "center": [ 0.1164, -0.2178 ] }, + { "ring": 1, "number": 27, "center": [ 0.1567, -0.1909 ] }, + { "ring": 1, "number": 28, "center": [ 0.1909, -0.1567 ] }, + { "ring": 1, "number": 29, "center": [ 0.2178, -0.1164 ] }, + { "ring": 1, "number": 30, "center": [ 0.2364, -0.0717 ] }, + { "ring": 1, "number": 31, "center": [ 0.2458, -0.0242 ] }, + { "ring": 2, "number": 0, "center": [ 0.2936, 0.0289 ] }, + { "ring": 2, "number": 1, "center": [ 0.2823, 0.0856 ] }, + { "ring": 2, "number": 2, "center": [ 0.2602, 0.1391 ] }, + { "ring": 2, "number": 3, "center": [ 0.2280, 0.1871 ] }, + { "ring": 2, "number": 4, "center": [ 0.1871, 0.2280 ] }, + { "ring": 2, "number": 5, "center": [ 0.1391, 0.2602 ] }, + { "ring": 2, "number": 6, "center": [ 0.0856, 0.2823 ] }, + { "ring": 2, "number": 7, "center": [ 0.0289, 0.2936 ] }, + { "ring": 2, "number": 8, "center": [ -0.0289, 0.2936 ] }, + { "ring": 2, "number": 9, "center": [ -0.0856, 0.2823 ] }, + { "ring": 2, "number": 10, "center": [ -0.1391, 0.2602 ] }, + { "ring": 2, "number": 11, "center": [ -0.1871, 0.2280 ] }, + { "ring": 2, "number": 12, "center": [ -0.2280, 0.1871 ] }, + { "ring": 2, "number": 13, "center": [ -0.2602, 0.1391 ] }, + { "ring": 2, "number": 14, "center": [ -0.2823, 0.0856 ] }, + { "ring": 2, "number": 15, "center": [ -0.2936, 0.0289 ] }, + { "ring": 2, "number": 16, "center": [ -0.2936, -0.0289 ] }, + { "ring": 2, "number": 17, "center": [ -0.2823, -0.0856 ] }, + { "ring": 2, "number": 18, "center": [ -0.2602, -0.1391 ] }, + { "ring": 2, "number": 19, "center": [ -0.2280, -0.1871 ] }, + { "ring": 2, "number": 20, "center": [ -0.1871, -0.2280 ] }, + { "ring": 2, "number": 21, "center": [ -0.1391, -0.2602 ] }, + { "ring": 2, "number": 22, "center": [ -0.0856, -0.2823 ] }, + { "ring": 2, "number": 23, "center": [ -0.0289, -0.2936 ] }, + { "ring": 2, "number": 24, "center": [ 0.0289, -0.2936 ] }, + { "ring": 2, "number": 25, "center": [ 0.0856, -0.2823 ] }, + { "ring": 2, "number": 26, "center": [ 0.1391, -0.2602 ] }, + { "ring": 2, "number": 27, "center": [ 0.1871, -0.2280 ] }, + { "ring": 2, "number": 28, "center": [ 0.2280, -0.1871 ] }, + { "ring": 2, "number": 29, "center": [ 0.2602, -0.1391 ] }, + { "ring": 2, "number": 30, "center": [ 0.2823, -0.0856 ] }, + { "ring": 2, "number": 31, "center": [ 0.2936, -0.0289 ] }, + { "ring": 3, "number": 0, "center": [ 0.3384, 0.0333 ] }, + { "ring": 3, "number": 1, "center": [ 0.3254, 0.0987 ] }, + { "ring": 3, "number": 2, "center": [ 0.2999, 0.1603 ] }, + { "ring": 3, "number": 3, "center": [ 0.2628, 0.2157 ] }, + { "ring": 3, "number": 4, "center": [ 0.2157, 0.2628 ] }, + { "ring": 3, "number": 5, "center": [ 0.1603, 0.2999 ] }, + { "ring": 3, "number": 6, "center": [ 0.0987, 0.3254 ] }, + { "ring": 3, "number": 7, "center": [ 0.0333, 0.3384 ] }, + { "ring": 3, "number": 8, "center": [ -0.0333, 0.3384 ] }, + { "ring": 3, "number": 9, "center": [ -0.0987, 0.3254 ] }, + { "ring": 3, "number": 10, "center": [ -0.1603, 0.2999 ] }, + { "ring": 3, "number": 11, "center": [ -0.2157, 0.2628 ] }, + { "ring": 3, "number": 12, "center": [ -0.2628, 0.2157 ] }, + { "ring": 3, "number": 13, "center": [ -0.2999, 0.1603 ] }, + { "ring": 3, "number": 14, "center": [ -0.3254, 0.0987 ] }, + { "ring": 3, "number": 15, "center": [ -0.3384, 0.0333 ] }, + { "ring": 3, "number": 16, "center": [ -0.3384, -0.0333 ] }, + { "ring": 3, "number": 17, "center": [ -0.3254, -0.0987 ] }, + { "ring": 3, "number": 18, "center": [ -0.2999, -0.1603 ] }, + { "ring": 3, "number": 19, "center": [ -0.2628, -0.2157 ] }, + { "ring": 3, "number": 20, "center": [ -0.2157, -0.2628 ] }, + { "ring": 3, "number": 21, "center": [ -0.1603, -0.2999 ] }, + { "ring": 3, "number": 22, "center": [ -0.0987, -0.3254 ] }, + { "ring": 3, "number": 23, "center": [ -0.0333, -0.3384 ] }, + { "ring": 3, "number": 24, "center": [ 0.0333, -0.3384 ] }, + { "ring": 3, "number": 25, "center": [ 0.0987, -0.3254 ] }, + { "ring": 3, "number": 26, "center": [ 0.1603, -0.2999 ] }, + { "ring": 3, "number": 27, "center": [ 0.2157, -0.2628 ] }, + { "ring": 3, "number": 28, "center": [ 0.2628, -0.2157 ] }, + { "ring": 3, "number": 29, "center": [ 0.2999, -0.1603 ] }, + { "ring": 3, "number": 30, "center": [ 0.3254, -0.0987 ] }, + { "ring": 3, "number": 31, "center": [ 0.3384, -0.0333 ] }, + { "ring": 4, "number": 0, "center": [ 0.3772, 0.0371 ] }, + { "ring": 4, "number": 1, "center": [ 0.3627, 0.1100 ] }, + { "ring": 4, "number": 2, "center": [ 0.3342, 0.1787 ] }, + { "ring": 4, "number": 3, "center": [ 0.2930, 0.2404 ] }, + { "ring": 4, "number": 4, "center": [ 0.2404, 0.2930 ] }, + { "ring": 4, "number": 5, "center": [ 0.1787, 0.3342 ] }, + { "ring": 4, "number": 6, "center": [ 0.1100, 0.3627 ] }, + { "ring": 4, "number": 7, "center": [ 0.0371, 0.3772 ] }, + { "ring": 4, "number": 8, "center": [ -0.0371, 0.3772 ] }, + { "ring": 4, "number": 9, "center": [ -0.1100, 0.3627 ] }, + { "ring": 4, "number": 10, "center": [ -0.1787, 0.3342 ] }, + { "ring": 4, "number": 11, "center": [ -0.2404, 0.2930 ] }, + { "ring": 4, "number": 12, "center": [ -0.2930, 0.2404 ] }, + { "ring": 4, "number": 13, "center": [ -0.3342, 0.1787 ] }, + { "ring": 4, "number": 14, "center": [ -0.3627, 0.1100 ] }, + { "ring": 4, "number": 15, "center": [ -0.3772, 0.0371 ] }, + { "ring": 4, "number": 16, "center": [ -0.3772, -0.0371 ] }, + { "ring": 4, "number": 17, "center": [ -0.3627, -0.1100 ] }, + { "ring": 4, "number": 18, "center": [ -0.3342, -0.1787 ] }, + { "ring": 4, "number": 19, "center": [ -0.2930, -0.2404 ] }, + { "ring": 4, "number": 20, "center": [ -0.2404, -0.2930 ] }, + { "ring": 4, "number": 21, "center": [ -0.1787, -0.3342 ] }, + { "ring": 4, "number": 22, "center": [ -0.1100, -0.3627 ] }, + { "ring": 4, "number": 23, "center": [ -0.0371, -0.3772 ] }, + { "ring": 4, "number": 24, "center": [ 0.0371, -0.3772 ] }, + { "ring": 4, "number": 25, "center": [ 0.1100, -0.3627 ] }, + { "ring": 4, "number": 26, "center": [ 0.1787, -0.3342 ] }, + { "ring": 4, "number": 27, "center": [ 0.2404, -0.2930 ] }, + { "ring": 4, "number": 28, "center": [ 0.2930, -0.2404 ] }, + { "ring": 4, "number": 29, "center": [ 0.3342, -0.1787 ] }, + { "ring": 4, "number": 30, "center": [ 0.3627, -0.1100 ] }, + { "ring": 4, "number": 31, "center": [ 0.3772, -0.0371 ] }, + { "ring": 5, "number": 0, "center": [ 0.4120, 0.0406 ] }, + { "ring": 5, "number": 1, "center": [ 0.3962, 0.1202 ] }, + { "ring": 5, "number": 2, "center": [ 0.3651, 0.1952 ] }, + { "ring": 5, "number": 3, "center": [ 0.3200, 0.2626 ] }, + { "ring": 5, "number": 4, "center": [ 0.2626, 0.3200 ] }, + { "ring": 5, "number": 5, "center": [ 0.1952, 0.3651 ] }, + { "ring": 5, "number": 6, "center": [ 0.1202, 0.3962 ] }, + { "ring": 5, "number": 7, "center": [ 0.0406, 0.4120 ] }, + { "ring": 5, "number": 8, "center": [ -0.0406, 0.4120 ] }, + { "ring": 5, "number": 9, "center": [ -0.1202, 0.3962 ] }, + { "ring": 5, "number": 10, "center": [ -0.1952, 0.3651 ] }, + { "ring": 5, "number": 11, "center": [ -0.2626, 0.3200 ] }, + { "ring": 5, "number": 12, "center": [ -0.3200, 0.2626 ] }, + { "ring": 5, "number": 13, "center": [ -0.3651, 0.1952 ] }, + { "ring": 5, "number": 14, "center": [ -0.3962, 0.1202 ] }, + { "ring": 5, "number": 15, "center": [ -0.4120, 0.0406 ] }, + { "ring": 5, "number": 16, "center": [ -0.4120, -0.0406 ] }, + { "ring": 5, "number": 17, "center": [ -0.3962, -0.1202 ] }, + { "ring": 5, "number": 18, "center": [ -0.3651, -0.1952 ] }, + { "ring": 5, "number": 19, "center": [ -0.3200, -0.2626 ] }, + { "ring": 5, "number": 20, "center": [ -0.2626, -0.3200 ] }, + { "ring": 5, "number": 21, "center": [ -0.1952, -0.3651 ] }, + { "ring": 5, "number": 22, "center": [ -0.1202, -0.3962 ] }, + { "ring": 5, "number": 23, "center": [ -0.0406, -0.4120 ] }, + { "ring": 5, "number": 24, "center": [ 0.0406, -0.4120 ] }, + { "ring": 5, "number": 25, "center": [ 0.1202, -0.3962 ] }, + { "ring": 5, "number": 26, "center": [ 0.1952, -0.3651 ] }, + { "ring": 5, "number": 27, "center": [ 0.2626, -0.3200 ] }, + { "ring": 5, "number": 28, "center": [ 0.3200, -0.2626 ] }, + { "ring": 5, "number": 29, "center": [ 0.3651, -0.1952 ] }, + { "ring": 5, "number": 30, "center": [ 0.3962, -0.1202 ] }, + { "ring": 5, "number": 31, "center": [ 0.4120, -0.0406 ] }, + { "ring": 6, "number": 0, "center": [ 0.4458, 0.0439 ] }, + { "ring": 6, "number": 1, "center": [ 0.4287, 0.1300 ] }, + { "ring": 6, "number": 2, "center": [ 0.3951, 0.2112 ] }, + { "ring": 6, "number": 3, "center": [ 0.3463, 0.2842 ] }, + { "ring": 6, "number": 4, "center": [ 0.2842, 0.3463 ] }, + { "ring": 6, "number": 5, "center": [ 0.2112, 0.3951 ] }, + { "ring": 6, "number": 6, "center": [ 0.1300, 0.4287 ] }, + { "ring": 6, "number": 7, "center": [ 0.0439, 0.4458 ] }, + { "ring": 6, "number": 8, "center": [ -0.0439, 0.4458 ] }, + { "ring": 6, "number": 9, "center": [ -0.1300, 0.4287 ] }, + { "ring": 6, "number": 10, "center": [ -0.2112, 0.3951 ] }, + { "ring": 6, "number": 11, "center": [ -0.2842, 0.3463 ] }, + { "ring": 6, "number": 12, "center": [ -0.3463, 0.2842 ] }, + { "ring": 6, "number": 13, "center": [ -0.3951, 0.2112 ] }, + { "ring": 6, "number": 14, "center": [ -0.4287, 0.1300 ] }, + { "ring": 6, "number": 15, "center": [ -0.4458, 0.0439 ] }, + { "ring": 6, "number": 16, "center": [ -0.4458, -0.0439 ] }, + { "ring": 6, "number": 17, "center": [ -0.4287, -0.1300 ] }, + { "ring": 6, "number": 18, "center": [ -0.3951, -0.2112 ] }, + { "ring": 6, "number": 19, "center": [ -0.3463, -0.2842 ] }, + { "ring": 6, "number": 20, "center": [ -0.2842, -0.3463 ] }, + { "ring": 6, "number": 21, "center": [ -0.2112, -0.3951 ] }, + { "ring": 6, "number": 22, "center": [ -0.1300, -0.4287 ] }, + { "ring": 6, "number": 23, "center": [ -0.0439, -0.4458 ] }, + { "ring": 6, "number": 24, "center": [ 0.0439, -0.4458 ] }, + { "ring": 6, "number": 25, "center": [ 0.1300, -0.4287 ] }, + { "ring": 6, "number": 26, "center": [ 0.2112, -0.3951 ] }, + { "ring": 6, "number": 27, "center": [ 0.2842, -0.3463 ] }, + { "ring": 6, "number": 28, "center": [ 0.3463, -0.2842 ] }, + { "ring": 6, "number": 29, "center": [ 0.3951, -0.2112 ] }, + { "ring": 6, "number": 30, "center": [ 0.4287, -0.1300 ] }, + { "ring": 6, "number": 31, "center": [ 0.4458, -0.0439 ] }, + { "ring": 7, "number": 0, "center": [ 0.4757, 0.0469 ] }, + { "ring": 7, "number": 1, "center": [ 0.4574, 0.1388 ] }, + { "ring": 7, "number": 2, "center": [ 0.4216, 0.2253 ] }, + { "ring": 7, "number": 3, "center": [ 0.3695, 0.3032 ] }, + { "ring": 7, "number": 4, "center": [ 0.3032, 0.3695 ] }, + { "ring": 7, "number": 5, "center": [ 0.2253, 0.4216 ] }, + { "ring": 7, "number": 6, "center": [ 0.1388, 0.4574 ] }, + { "ring": 7, "number": 7, "center": [ 0.0469, 0.4757 ] }, + { "ring": 7, "number": 8, "center": [ -0.0469, 0.4757 ] }, + { "ring": 7, "number": 9, "center": [ -0.1388, 0.4574 ] }, + { "ring": 7, "number": 10, "center": [ -0.2253, 0.4216 ] }, + { "ring": 7, "number": 11, "center": [ -0.3032, 0.3695 ] }, + { "ring": 7, "number": 12, "center": [ -0.3695, 0.3032 ] }, + { "ring": 7, "number": 13, "center": [ -0.4216, 0.2253 ] }, + { "ring": 7, "number": 14, "center": [ -0.4574, 0.1388 ] }, + { "ring": 7, "number": 15, "center": [ -0.4757, 0.0469 ] }, + { "ring": 7, "number": 16, "center": [ -0.4757, -0.0469 ] }, + { "ring": 7, "number": 17, "center": [ -0.4574, -0.1388 ] }, + { "ring": 7, "number": 18, "center": [ -0.4216, -0.2253 ] }, + { "ring": 7, "number": 19, "center": [ -0.3695, -0.3032 ] }, + { "ring": 7, "number": 20, "center": [ -0.3032, -0.3695 ] }, + { "ring": 7, "number": 21, "center": [ -0.2253, -0.4216 ] }, + { "ring": 7, "number": 22, "center": [ -0.1388, -0.4574 ] }, + { "ring": 7, "number": 23, "center": [ -0.0469, -0.4757 ] }, + { "ring": 7, "number": 24, "center": [ 0.0469, -0.4757 ] }, + { "ring": 7, "number": 25, "center": [ 0.1388, -0.4574 ] }, + { "ring": 7, "number": 26, "center": [ 0.2253, -0.4216 ] }, + { "ring": 7, "number": 27, "center": [ 0.3032, -0.3695 ] }, + { "ring": 7, "number": 28, "center": [ 0.3695, -0.3032 ] }, + { "ring": 7, "number": 29, "center": [ 0.4216, -0.2253 ] }, + { "ring": 7, "number": 30, "center": [ 0.4574, -0.1388 ] }, + { "ring": 7, "number": 31, "center": [ 0.4757, -0.0469 ] }, + { "mark": 1, "number": 0, "center": [ 0.3141, -0.1301 ] }, + { "mark": 1, "number": 1, "center": [ 0.0933, -0.4688 ] }, + ] +} diff --git a/mask/json.c b/mask/json.c new file mode 100644 index 0000000..aa251b0 --- /dev/null +++ b/mask/json.c @@ -0,0 +1,205 @@ +#include +#include + +#define FREE(ptr) do{free(ptr); ptr = NULL;}while(0) + +typedef struct{ + float x0; // left border + float y0; // lower border + float w; // width + float h; // height +} BBox; +typedef enum{ + H_SQUARE // square hole + ,H_ELLIPSE // elliptic hole + ,H_UNDEF +} HoleType; +typedef struct{ + BBox box; // bounding box of hole + int type; // type, in case of round hole borders of box are tangents to hole +} aHole; + +aHole globHole; + +double get_jdouble(json_object *jobj){ + enum json_type type = json_object_get_type(jobj); + double val; + switch(type){ + case json_type_double: + val = json_object_get_double(jobj); + break; + case json_type_int: + val = json_object_get_int(jobj); + break; + default: + fprintf(stderr, "Wrong value! Get non-number!\n"); + exit(-1); + } + return val; +} + +double *json_get_array(json_object *jobj, int *getLen){ + enum json_type type; + json_object *jarray = jobj; + int arraylen = json_object_array_length(jarray); + double *arr = calloc(arraylen, sizeof(double)); + int L = 0; + int i; + json_object *jvalue; + for (i=0; i< arraylen; i++){ + jvalue = json_object_array_get_idx(jarray, i); + type = json_object_get_type(jvalue); + if(type == json_type_array){ // nested arrays is error + fprintf(stderr, "Invalid file format! Found nested arrays!\n"); + exit(-1); + } + else if (type != json_type_object){ + arr[L++] = get_jdouble(jvalue); + } + else{ // non-numerical data? + fprintf(stderr, "Invalid file format! Non-numerical data in array!\n"); + exit(-1); + } + } + if(L == 0) FREE(arr); + if(getLen) *getLen = L; + return arr; +} + +void get_obj_params(json_object *jobj, aHole *H){ + double *arr = NULL; + int Len; + enum json_type type; + if(!H){ + fprintf(stderr, "Error: NULL instead of aHole structure!\n"); + exit(-1); + } + memcpy(H, &globHole, sizeof(aHole)); // initialize hole by global values + json_object *o = json_object_object_get(jobj, "shape"); + if(o){ + const char *ptr = json_object_get_string(o); + if(strcmp(ptr, "square") == 0) H->type = H_SQUARE; + else if(strcmp(ptr, "round") == 0 || strcmp(ptr, "ellipse") == 0 ) H->type = H_ELLIPSE; + else H->type = H_UNDEF; + } + o = json_object_object_get(jobj, "radius"); + if(o){ + type = json_object_get_type(o); + if(type == json_type_int || type == json_type_double){ // circle / square + double R = json_object_get_double(o); + H->box.w = H->box.h = R * 2.; + }else if(type == json_type_array){ // ellipse / rectangle + if(!(arr = json_get_array(o, &Len)) || Len != 2){ + fprintf(stderr, "\"radius\" array must consist of two doubles!\n"); + exit(-1); + } + H->box.w = arr[0] * 2.; H->box.h = arr[1] * 2.; + FREE(arr); + }else{ + fprintf(stderr, "\"radius\" must be a number or an array of two doubles!\n"); + exit(-1); + } + } + o = json_object_object_get(jobj, "center"); + if(o){ + if(!(arr = json_get_array(o, &Len)) || Len != 2){ + fprintf(stderr, "\"center\" must contain an array of two doubles!\n"); + exit(-1); + } + H->box.x0 = arr[0] - H->box.w/2.; + H->box.y0 = arr[1] - H->box.h/2.; + FREE(arr); + }else{ + o = json_object_object_get(jobj, "bbox"); + if(o){ + if(!(arr = json_get_array(o, &Len)) || Len != 4){ + fprintf(stderr, "\"bbox\" must contain an array of four doubles!\n"); + exit(-1); + } + H->box.x0 = arr[0]; H->box.y0 = arr[1]; H->box.w = arr[2]; H->box.h = arr[3]; + FREE(arr); + } + } +} + +aHole *json_parse_holesarray(json_object *jobj, int *getLen){ + enum json_type type; + json_object *jarray = jobj; + int arraylen = json_object_array_length(jarray), i; + aHole *H = calloc(arraylen, sizeof(aHole)); + json_object *jvalue; + for (i=0; i < arraylen; i++){ + jvalue = json_object_array_get_idx(jarray, i); + type = json_object_get_type(jvalue); + if(type == json_type_object){ + get_obj_params(jvalue, &H[i]); + }else{ + fprintf(stderr, "Invalid holes array format!\n"); + exit(-1); + } + } + if(getLen) *getLen = arraylen; + return H; +} + +char *gettype(aHole *H){ + char *ret; + switch(H->type){ + case H_SQUARE: + ret = "square"; + break; + case H_ELLIPSE: + ret = "ellipse"; + break; + default: + ret = "undefined"; + } + return ret; +} + +int main(int argc, char **argv){ + char *ptr; + struct stat statbuf; + enum json_type type; + int fd, i, HolesNum; + aHole *HolesArray; + size_t Mlen; + if(argc == 2){ + if ((fd = open (argv[1], O_RDONLY)) < 0) err (1, "Can't open %s for reading", argv[1]); + if (fstat (fd, &statbuf) < 0) err (1, "Fstat error"); + Mlen = statbuf.st_size; + if ((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED) + err(1, "Mmap error for input"); + }else{ + fprintf(stderr, "No file given!\n"); + exit(-1); + } + json_object * jobj = json_tokener_parse(ptr); + get_obj_params(jobj, &globHole); + json_object *o = json_object_object_get(jobj, "holes"); + if(!o){ + fprintf(stderr, "Corrupted file: no holes found!\n"); + exit(-1); + } + type = json_object_get_type(o); + if(type == json_type_object){ // single hole + HolesArray = calloc(1, sizeof(aHole)); + assert(HolesArray); + HolesNum = 1; + get_obj_params(o, HolesArray); + }else{ // array of holes + HolesArray = json_parse_holesarray(o, &HolesNum); + } + if(!HolesArray || HolesNum < 1){ + fprintf(stderr, "Didn't find any holes in json file!\n"); + exit(-1); + } + printf("Readed %d holes\n", HolesNum); + for(i = 0; i < HolesNum; i++){ + BBox *B = &HolesArray[i].box; + printf("Hole %d: type=%s, bbox={%g, %g, %g, %g}\n", i, gettype(&HolesArray[i]), + B->x0, B->y0, B->w, B->h); + } + munmap(ptr, Mlen); + return 0; +} diff --git a/mask/make_mask.m b/mask/make_mask.m new file mode 100644 index 0000000..a164145 --- /dev/null +++ b/mask/make_mask.m @@ -0,0 +1,34 @@ + +function make_mask() +% построение гартмановской маски +% SS - размер маски + f = fopen("holes.json", "w"); + R = [175 247 295 340 379 414 448 478] * 1e-3; % радиусы колец на гартманограмме + HoleR = 7.5e-3; % радиус отверстий - 7.5мм + R0 = .6; % радиус самой гартманограммы + alpha0 = pi/32; % смещение лучей относительно секущих горизонтали/вертикали + Angles = [0:31] * 2 * alpha0 + alpha0; % углы, по которым располагаются лучи + % для того, чтобы разместить на маске окружности, создадим маску + % окружности: zeros(15) с единицами там, где должна быть дырка. Затем + % пометим единицами в mask те точки, куда должен попадать левый верхний + fprintf(f, "{\n\t\"shape\": \"round\", \"radius\": %f,\n\t\"holes\": [\n" , HoleR); + for i = [1 : size(R,2)] % цикл по кольцам + x = R(i) * cos(Angles); + y = R(i) * sin(Angles); + %fprintf(f, "\t\t{\"ring\": %d, \"center\": [%f, %f]\n", i, x[j], y[j]]); + printR(f, sprintf("\"ring\": %d", i-1), x, y); + endfor + % помечаем маркеры + x = R([4 8]) .* cos([-2 -7]* 2 * alpha0); + y = R([4 8]) .* sin([-2 -7]* 2 * alpha0); + %fprintf(f, "\t\t{\"marker\", \"center\": [%f, %f]\n", x, y); + printR(f, sprintf("\"mark\": 1"), x, y); + fprintf(f, "\t]\n}\n"); + fclose(f); +endfunction + +function printR(f, msg, x, y) + for i = 1:size(x,2) + fprintf(f, "\t\t{ %9s, \"number\": %2d, \"center\": [ %7.4f, %7.4f ] },\n", msg, i-1, x(i), y(i)); + endfor +endfunction diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..c406e20 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,107 @@ +cmake_minimum_required(VERSION 2.6) +if(DEFINED DEBUG) + add_definitions(-DEBUG) +endif() +aux_source_directory(. SOURCES) +set(CUFILE CUDA.cu) +set(CFLAGS -O3 -Wall -Werror -std=gnu99) +set(LCPATH ${CMAKE_CURRENT_SOURCE_DIR}/locale/ru) +set(PO_FILE ${LCPATH}/messages.po) +set(MO_FILE ${LCPATH}/LC_MESSAGES/${PROJ}.mo) +set(RU_FILE ${LCPATH}/ru.po) +find_package(PkgConfig REQUIRED) +#include(FindOpenMP) +find_package(OpenMP) +find_package(JPEG) +find_package(TIFF) +find_package(PNG) +#find_package(OpenGL REQUIRED) +#find_package(GTK2 REQUIRED) +pkg_check_modules(${PROJ} REQUIRED +# gtkglext-1.0>=0.7.0 +# gtkglext-x11-1.0>=0.7.0 + cfitsio>=3.0 + json>=0.9 +# libpng>=1.5 +# fftw3>=3.2.0 + ) +if(OPENMP_FOUND) + add_definitions(-DOMP) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") +endif() +list(APPEND ${PROJ}_INCLUDE_DIRS ${JPEG_INCLUDE_DIR} ${PNG_INCLUDE_DIR} ${TIFF_INCLUDE_DIR}) +list(APPEND ${PROJ}_LIBRARIES ${JPEG_LIBRARY} ${PNG_LIBRARY} ${TIFF_LIBRARY}) +add_definitions(-D_XOPEN_SOURCE=1000 -D__JPEG=${JPEG_FOUND} -D__TIFF=${TIFF_FOUND} -D__PNG=${PNG_FOUND}) +#list(APPEND ${${PROJ}_LIBRARY_DIRS} ) +#set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -lfftw3_threads") +#if(NOT DEFINED NO_LEPTONICA) +# pkg_check_modules(LIBLEPT liblept) +#endif() +#if(NOT DEFINED NO_GSL) +# pkg_check_modules(GSL gsl) +#endif() +#if(NOT DEFINED GSL_VERSION) +# message("GSL not found, some mathematics functions wouldn't be avialable") +#else() +# add_definitions(-DGSL_FOUND) +# list(APPEND ${${PROJ}_LIBRARIES} ${GSL_LIBRARIES}) +# list(APPEND ${${PROJ}_INCLUDE_DIRS} ${GSL_INCLUDE_DIRS}) +# list(APPEND ${${PROJ}_LIBRARY_DIRS} ${GSL_LIBRARY_DIRS}) +#endif() +#if(NOT DEFINED LIBLEPT_VERSION) +# message("Leptonica library not found, some functions wouldn't be avialable") +#else() +# add_definitions(-DLEPTONICA_FOUND) +# list(APPEND ${${PROJ}_LIBRARIES} ${LIBLEPT_LIBRARIES}) +# list(APPEND ${${PROJ}_INCLUDE_DIRS} ${LIBLEPT_INCLUDE_DIRS}) +# list(APPEND ${${PROJ}_LIBRARY_DIRS} ${LIBLEPT_LIBRARY_DIRS}) +#endif() + +include_directories(include ${${PROJ}_INCLUDE_DIRS}) +link_directories(${${PROJ}_LIBRARY_DIRS}) +add_definitions(-DPACKAGE_VERSION=\"${PROJ_VER}\" -DGETTEXT_PACKAGE=\"${PROJ}\" + -DLOCALEDIR=\"${LOCALEDIR}\" ${CFLAGS}) + +if(CUDA_FOUND) + set(TEST ${CMAKE_BINARY_DIR}/test) + set(TESTSRC ${CMAKE_CURRENT_SOURCE_DIR}/test/capability.cu) + execute_process(COMMAND nvcc -lcuda ${TESTSRC} -o ${TEST}) + execute_process(COMMAND ${TEST} OUTPUT_VARIABLE CUDA_ARCH) + message("Cuda architecture: ${CUDA_ARCH}") + list(APPEND CUDA_NVCC_FLAGS --use_fast_math ${CUDA_ARCH}) + #-gencode arch=compute_20,code=sm_20 -gencode arch=compute_12,code=sm_12 -gencode arch=compute_13,code=sm_13) + cuda_include_directories(include) + cuda_add_executable(../${PROJ} ${SOURCES} ${CUFILE} ${PO_FILE} ${MO_FILE}) + target_link_libraries( ../${PROJ} ${${PROJ}_LIBRARIES} ${CUDA_curand_LIBRARY} + #${CUDA_CUFFT_LIBRARIES} + ) +else(CUDA_FOUND) + add_executable(../${PROJ} ${SOURCES} ${PO_FILE} ${MO_FILE}) + target_link_libraries( ../${PROJ} ${${PROJ}_LIBRARIES} + ${CMAKE_THREAD_LIBS_INIT} + ) +endif(CUDA_FOUND) + +# Installation of the program +INSTALL(FILES ${MO_FILE} DESTINATION ${CMAKE_INSTALL_LOCALEDIR}) +INSTALL(TARGETS ../${PROJ} DESTINATION "bin") + +#if(DEFINED DEBUG) + 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() +add_custom_command( + OUTPUT ${PO_FILE} + COMMAND ${GETTEXT_XGETTEXT_EXECUTABLE} -D ${CMAKE_CURRENT_SOURCE_DIR} --from-code=koi8-r ${SOURCES} -c -k_ -kN_ -o ${PO_FILE} + COMMAND sed 's/charset=UTF-8/charset=koi8-r/' ${PO_FILE} | enconv > tmp && mv -f tmp ${PO_FILE} + COMMAND ${GETTEXT_MSGMERGE_EXECUTABLE} -Uis ${RU_FILE} ${PO_FILE} + DEPENDS ${SOURCES}) +add_custom_command( + OUTPUT ${MO_FILE} + COMMAND ${GETTEXT_MSGFMT_EXECUTABLE} ${RU_FILE} -o ${MO_FILE} + DEPENDS ${RU_FILE}) +#endif(DEFINED DEBUG) diff --git a/src/CPU.c b/src/CPU.c new file mode 100644 index 0000000..9318fb9 --- /dev/null +++ b/src/CPU.c @@ -0,0 +1,113 @@ +/* + * CPU.c - CPU variants (if possible - on OpenMP) of main functions + * + * 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. + */ + +#define CPU_C +#include "mkHartmann.h" +#include "wrapper.h" + +// RANDOM NUMBER generation block ============================================> +static int rand_initialized = 0; +/* + * Fills linear array of size sz by random float numbers + * name: fillrandarr + * @param sz - size of array + * @param arr - array data (allocated outside) + * @param Amp - amplitude of array + * @return 0 if failed + */ +int CPUfillrandarr(size_t sz, float *arr, float Amp){ + FNAME(); + if(!rand_initialized){ + srand48(throw_random_seed()); + rand_initialized = 1; + } + size_t i; + OMP_FOR() + for(i = 0; i < sz; i++) + arr[i] = (float)drand48() * Amp; + return 1; +} +// <===================================== end of random number generation block + +// BICUBIC interpolation block ===============================================> +inline float p_ta(float t, float t2, float t3, + float a,float b,float c,float d){ + return (2*b + t*(-a+c) + t2*(2*a-5*b+4*c-d) + t3*(-a+3*b-3*c+d)) / 2.f; +} + +int CPUbicubic_interp(float *out, float *in, + size_t oH, size_t oW, size_t iH, size_t iW){ + FNAME(); + float fracX = (float)(iW-1)/(oW-1), fracY = (float)(iH-1)/(oH-1); + size_t X, Y, ym1,y1,y2, xm1,x1,x2; // pixel coordinates on output + size_t Ym = iH - 1, Xm = iW - 1, Pcur; + float x,y; // coordinates on output in value of input + OMP_FOR() + for(Y = 0; Y < oH; Y++){ + // we can't do "y+=fracY" because of possible cumulative error + y = (float)Y * fracY; + int y0 = floor(y); + float pty = y - (float)y0, pty2 = pty*pty, pty3 = pty*pty2; + ym1 = y0-1; y1 = y0 + 1; y2 = y0 + 2; + if(y0 == 0) ym1 = 0; + if(y1 > Ym) y1 = Ym; + if(y2 > Ym) y2 = Ym; + for(X = 0, Pcur = Y * oW; X < oW; X++, Pcur++){ + // we can't do "x+=fracX" because of possible cumulative error + x = (float)X * fracX; + int x0 = floor(x); + float ptx = x - (float)x0, ptx2 = ptx*ptx, ptx3 = ptx*ptx2; + xm1 = x0-1; x1 = x0 + 1; x2 = x0 + 2; + if(x0 == 0) xm1 = 0; + if(x1 > Xm) x1 = Xm; + if(x2 > Xm) x2 = Xm; + #define TEX(x,y) (in[iW*y + x]) + #define TX(y) p_ta(ptx, ptx2, ptx3, TEX(xm1,y), \ + TEX(x0,y), TEX(x1,y), TEX(x2,y)) + out[Pcur] = p_ta( + pty, pty2, pty3, + TX(ym1), TX(y0), TX(y1), TX(y2)); + #undef TX + #undef TEX + } + } + return 1; +} +// <========================================= end of BICUBIC interpolation block + +int CPUmkmirDeviat(float *map _U_, size_t mapWH _U_, float mirSZ _U_, + mirDeviations * mirDev _U_){ + FNAME(); + return 0; +} + +int CPUgetPhotonXY(float *xout _U_, float *yout _U_, int R _U_ ,mirDeviations *D _U_, + mirPar *mirParms _U_, size_t N_photons _U_, BBox *box _U_){ + FNAME(); + return 0; +} + +/* +int CPUfillImage(float *phX _U_, float *phY _U_, size_t ph_sz _U_, + float *image _U_, size_t imW _U_, size_t imH _U_, BBox *imbox _U_){ + return 0; +} +*/ diff --git a/src/CUDA.cu b/src/CUDA.cu new file mode 100644 index 0000000..3b46d6e --- /dev/null +++ b/src/CUDA.cu @@ -0,0 +1,807 @@ +/* + * CUDA.cu - subroutines for GPU + * + * 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 "cutil_math.h" + + +//#include +//#include +//#include + +#define EXTERN extern "C" +#define CUDA_CU +#include "mkHartmann.h" +#include "wrapper.h" + +int SHMEMSZ = 16383; // default constants, changed runtime +int QBLKSZ = 16; // QBLKSZ = sqrt(LBLKSZ) +int LBLKSZ = 512; + +cudaError_t CUerr; +inline int CUERROR(char *str){ + if(CUerr != cudaSuccess){ + WARN("%s, %s", str, cudaGetErrorString(CUerr)); + return 1; + }else return 0; +} + +// Default error macro (return fail) +#define RETMACRO do{return 0;} while(0) +/* + * memory macros + * + * Each function using them must return an int type: + * 1 - in case of success + * 0 - in case of fail + * + * Function that use them should check this status and if fails, + * call CPU-based function for same computations + */ +static int ret; +/* + * this macro is for functions with alloc + * All functions should in their beginning call "ret = 1;" + * and in the end they should have a label "free_all:" after which + * located memory free operators + * in the end of functions should be "return ret;" + */ +#define FREERETMACRO do{ret = 0; goto free_all;} while(0) + +#define CUALLOC(var, size) do{ \ + CUerr = cudaMalloc((void**)&var, size); \ + if(CUERROR("CUDA: can't allocate memory")){ \ + FREERETMACRO; \ +}}while(0) +#define CUALLOCPITCH(var, p, W, H) do{ \ + CUerr = cudaMallocPitch((void**)&var, p, W, H);\ + if(CUERROR("CUDA: can't allocate memory")){ \ + FREERETMACRO; \ +}}while(0) +#define CUMOV2DEV(dest, src, size) do{ \ + CUerr = cudaMemcpy(dest, src, size, \ + cudaMemcpyHostToDevice); \ + if(CUERROR("CUDA: can't copy data to device")){\ + FREERETMACRO; \ +}}while(0) +#define CUMOV2DEVPITCH(dst, dp, src, w, h) do{ \ + CUerr = cudaMemcpy2D(dst, dp, src, w, w, h,\ + cudaMemcpyHostToDevice); \ + if(CUERROR("CUDA: can't copy data to device")){\ + FREERETMACRO; \ +}}while(0) +#define CUMOV2HOST(dest, src, size) do{ \ + CUerr = cudaMemcpy(dest, src, size, \ + cudaMemcpyDeviceToHost); \ + if(CUERROR("CUDA: can't copy data to host")){\ + FREERETMACRO; \ +}}while(0) +#define CUMOV2HOSTPITCH(dst,src,spitch,w,h) do{ \ + CUerr = cudaMemcpy2D(dst,w,src,spitch,w,h, \ + cudaMemcpyDeviceToHost); \ + if(CUERROR("CUDA: can't copy data to device")){\ + FREERETMACRO; \ +}}while(0) +#define CUFREE(var) do{cudaFree(var); var = NULL; }while(0) +#define CUFFTCALL(fn) do{ \ + cufftResult fres = fn; \ + if(CUFFT_SUCCESS != fres){ \ + WARN("CUDA fft error %d", fres); \ + FREERETMACRO; \ +}}while(0) + +texture cuTex;//1, cuTex2, cuTex3; +#define CUTEXTURE(t, data, W, H, pitch) do{ \ + CUerr = cudaBindTexture2D(NULL, t, \ + data, W, H, pitch); \ + if(CUERROR("CUDA: can't bind texture")){ \ + FREERETMACRO;}else{ \ + t.addressMode[0] = cudaAddressModeClamp; \ + t.addressMode[1] = cudaAddressModeClamp; \ + t.normalized = false; \ + t.filterMode = cudaFilterModePoint; \ +}}while(0) + +//#define _TEXTURE_(N, ...) CUTEXTURE(cuTex ## N, __VA_ARGS__) +//#define TEXDATA(...) _TEXTURE_(__VA_ARGS__) +//#define TEXTURE(N) cuTex ## N +#define TEXDATA(...) CUTEXTURE(cuTex, __VA_ARGS__) +#define TEXTURE() cuTex + +/** + * getting the videocard parameters + * @return 0 if check failed + */ +EXTERN int CUgetprops(){ + cudaDeviceProp dP; + CUdevice dev; + CUcontext ctx; + if(cudaSuccess != cudaGetDeviceProperties(&dP, 0)) return 0; + if(CUDA_SUCCESS != cuDeviceGet(&dev,0)) return 0; + // create context for program run: + if(CUDA_SUCCESS != cuCtxCreate(&ctx, 0, dev)) return 0; + printf("\nDevice: %s, totalMem=%zd, memPerBlk=%zd,\n", dP.name, dP.totalGlobalMem, dP.sharedMemPerBlock); + printf("warpSZ=%d, TPB=%d, TBDim=%dx%dx%d\n", dP.warpSize, dP.maxThreadsPerBlock, + dP.maxThreadsDim[0],dP.maxThreadsDim[1],dP.maxThreadsDim[2]); + printf("GridSz=%dx%dx%d, MemovrLap=%d, GPUs=%d\n", dP.maxGridSize[0], + dP.maxGridSize[1],dP.maxGridSize[2], + dP.deviceOverlap, dP.multiProcessorCount); + printf("canMAPhostMEM=%d\n", dP.canMapHostMemory); + printf("compute capability"); + green(" %d.%d.\n\n", dP.major, dP.minor); + if(dP.major > 1){ + SHMEMSZ = 49151; QBLKSZ = 32; LBLKSZ = 1024; + } + // cuCtxDetach(ctx); + return 1; +} + +/** + * check whether there is enough memory + * @param memsz - memory needed + * @param Free - free memory (maybe NULL if not needed) + * @param Total - total memory (maybe NULL if not needed) + * @return 0 if check failed + */ +EXTERN int CUgetMEM(size_t memsz, size_t *Free, size_t *Total){ + size_t theFree = 0, theTotal = 0; + if(CUDA_SUCCESS != cuMemGetInfo( &theFree, &theTotal )) return 0; + if(Free) *Free = theFree; + if(Total) *Total = theTotal; + if(theFree < memsz) return 0; + return 1; +} + +/* +size_t theFree, theTotal; +CUgetMEM(0, &theFree, &theTotal); +printf(_("MEMORY: free = ")); +green("%zdMB,", theFree / MB); +printf(_(" total= ")); +green("%zdMB\n", theTotal / MB); +*/ + +/** + * Memory allocation & initialisation test + * @param memsz - wanted memory to allocate + * @return 0 if check failed + */ +EXTERN int CUallocaTest(size_t memsz){ + char *mem; + ret = 1; + printf("cudaMalloc(char, %zd)\n", memsz); + CUALLOC(mem, memsz); + if(cudaSuccess != cudaMemset(mem, 0xaa, memsz)) ret = 0; +free_all: + CUFREE(mem); + return ret; +} + +/* + * + * + * RANDOM NUMBER generation block ============================================> + * + * + */ +/* + * Possible functions: + * curandGenerateUniform(curandGenerator_t generator, float *outputPtr, size_t num) + * curandGenerateNormal(curandGenerator_t generator, float *outputPtr, size_t n, float mean, float stddev) + * curandGenerateLogNormal(curandGenerator_t generator, float *outputPtr, size_t n, float mean, float stddev) + * curandGeneratePoisson(curandGenerator_t generator, unsigned int *outputPtr, size_t n, double lambda) + */ +#define CURAND_CALL(x) do { \ + curandStatus_t st = x; \ + if(st!=CURAND_STATUS_SUCCESS) { \ + WARN("CURAND error %d", st); \ + FREERETMACRO; \ +}} while(0) + + +__global__ void multiply_rand(float *arr, size_t sz, float Amp){ + size_t IDX = threadIdx.x + blockDim.x * blockIdx.x; + if(IDX >= sz) return; + arr[IDX] *= Amp; +} +static int rand_initialized = 0; +static curandGenerator_t gen; +/** + * Fills linear array of size sz by random float numbers + * @param sz - size of array + * @param arr - array data (allocated outside) + * @param Amp - amplitude of array + * @return 0 if failed + */ +EXTERN int CUfillrandarr(size_t sz, float *arr, float Amp){ + FNAME(); + ret = 1; + float *devmem; + if(!rand_initialized){ + CURAND_CALL(curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT)); + CURAND_CALL(curandSetPseudoRandomGeneratorSeed(gen, throw_random_seed())); + rand_initialized = 1; + } + CUALLOC(devmem, sz*sizeof(float)); + CURAND_CALL(curandGenerateUniform(gen, devmem, sz)); + cudaThreadSynchronize(); + if(fabs(Amp < 1.f) > FLT_EPSILON){ + size_t dimens = (sz + LBLKSZ - 1) / LBLKSZ; + multiply_rand<<>>(devmem, sz, Amp); + cudaThreadSynchronize(); + } + CUMOV2HOST(arr, devmem, sz*sizeof(float)); +free_all: + //CURAND_CALL(curandDestroyGenerator(gen)); + CUFREE(devmem); + return ret; +} +/* + * + * + * <===================================== end of random number generation block + * + * + */ + +/* + * + * + * BICUBIC interpolation block ===============================================> + * + * + */ +/** + * matrix operation: + * + * p(t, a, b, c, d) = + * / 0 2 0 0 \ / a \ + * = 1/2 ( 1 t t^2 t^3 ) * | -1 0 1 0 | * | b | + * | 2 -5 4 -1 | | c | + * \ -1 3 -3 1 / \ d / + */ +inline __device__ float p_ta(float t, float t2, float t3, + float a,float b,float c,float d){ + return (2*b + t*(-a+c) + t2*(2*a-5*b+4*c-d) + t3*(-a+3*b-3*c+d)) / 2.f; +} +// bicubic interpolation of texture in point with coordinates (x, y) +__device__ float interpolate_bicubic(float x, float y){ + #define TEX(X,Y) tex2D(TEXTURE(), X, y0+(Y)) + #define TX(Y) p_ta(pt.x, pt2.x, pt3.x, TEX(x0-1,Y), \ + TEX(x0,Y), TEX(x0+1,Y), TEX(x0+2,Y)) + float2 crd = make_float2(x,y); + float2 zeropt = floor(crd); + float2 pt = crd - zeropt; + float2 pt2 = pt*pt; + float2 pt3 = pt*pt2; + float x0 = zeropt.x, y0 = zeropt.y; + //return x; + return p_ta( + pt.y, pt2.y, pt3.y, + TX(-1), TX(0), TX(1), TX(2) + ); + #undef TEX + #undef TX +} + +/** + * calculate base coordinates for new image in point (X,Y) relative to old image + * @param out - new image data + * @param fracX, fracY - scale of new image pixel relative to old image + * @parami oW, oH - new image dimensions + */ +__global__ void bicubic_interp(float *out, + float fracX, float fracY, + unsigned int oW, unsigned int oH){ + int X, Y; // pixel coordinates on output + float x,y; // coordinates on output in value of input + X = threadIdx.x + blockDim.x * blockIdx.x; + Y = threadIdx.y + blockDim.y * blockIdx.y; + if(X >= oW || Y >= oH) return; + x = (float)X * fracX; + y = (float)Y * fracY; + out[Y*oW+X] = interpolate_bicubic(x, y); +} + +/** + * Interpolation of image by bicubic splines + * @param out_dev - new image data in device memory + * @param in - old image data + * @param oH, oW - new image size + * @param iH, iW - old image size + * @return 0 if fails + */ +EXTERN int CUbicubic_interpDEV(float **out_dev, float *in, + size_t oH, size_t oW, size_t iH, size_t iW){ + FNAME(); + size_t iWp = iW*sizeof(float), oSz = oH*oW*sizeof(float), pitch; + float *in_dev = NULL; + ret = 1; + dim3 blkdim(QBLKSZ, QBLKSZ); + dim3 griddim((oW+QBLKSZ-1)/QBLKSZ, (oH+QBLKSZ-1)/QBLKSZ); + CUALLOCPITCH(in_dev, &pitch, iWp, iH); + CUMOV2DEVPITCH(in_dev, pitch, in, iWp, iH); + CUALLOC(*out_dev, oSz); + TEXDATA(in_dev, iW, iH, pitch); + bicubic_interp<<>>(*out_dev, + (float)(iW-1)/(oW-1), (float)(iH-1)/(oH-1), oW, oH); + cudaThreadSynchronize(); +free_all: + cudaUnbindTexture(TEXTURE()); + CUFREE(in_dev); + return ret; +} +/** + * Interpolation of image by bicubic splines + * @param out, in - new and old images data + * @param oH, oW - new image size + * @param iH, iW - old image size + * @return 0 if fails + */ +EXTERN int CUbicubic_interp(float *out, float *in, + size_t oH, size_t oW, size_t iH, size_t iW){ + FNAME(); + float *out_dev = NULL; + ret = 1; + if(!CUbicubic_interpDEV(&out_dev,in,oH,oW,iH,iW)){FREERETMACRO;} + CUMOV2HOST(out, out_dev, oH*oW*sizeof(float)); +free_all: + CUFREE(out_dev); + return ret; +} +/* + * + * + * <========================================= end of BICUBIC interpolation block + * + * + */ + + +/* + * + * + * Begin of mirror normales block =============================================> + * + * + */ + + +/** + * Compute X & Y components of "curved mirror" gradient + * This components (interpolated) should be added to SURFACE gradient **before** normalisation + * @param Sz - image size in both directions + * @param mirPixSZ - pixel size in meters + * @param dZdX, dZdY - addition gradient components + * + * + * Normale of "curved mirror" computes so: + * N1 = -i*x/2f -j*y/2f + k ==> direction vector of ideal mirror + * dN ==> bicubic interpolation of gradient matrix of mirror surface deviations + * dN = i*dX + j*dY + * N = N1 + dN ==> "real" direction vector + * |N| = length(N) ==> its length + * n = N / |N| ==> normale to "real" mirror surface + * the normale n can be transformed due to mirror inclination and displation + * + * Let f be a direction vector of falling photon, then a=[-2*dot(f,n)*n] + * will be a twice projection of f to n with "-" sign (dot(x,y) == scalar product) + * reflected direction vector will be + * r = a + f ==> + * r = f - 2*dot(f,n)*n + * r doesn't need normalisation + * + * If (x0,y0,z0) is mirror surface coordinates, z0 = z(x,y) + dZ, wehere dZ interpolated; + * then image coordinates on plane z=Zx would be: + * X = x0 + k*rx, Y = y0 + k*ry, k = (Zx - z0) / rz. + * + * First check point on a mask plane; if photon falls to a hole in mask, increment + * value on corresponding image pixel + */ +// 1./sqrt(2) +#define DIVSQ2 0.7071068f +// 2*(1+2/sqrt(2)) +#define WEIGHT 482.842712f +__global__ void calcMir_dXdY(size_t Sz, float mirPixSZ, + float *dZdX, float *dZdY){ + #define TEX(X,Y) tex2D(TEXTURE(), X0+(X), Y0+(Y))*100.f + #define PAIR(X,Y) (TEX((X),(Y))-TEX(-(X),-(Y))) + int X0,Y0; + X0 = threadIdx.x + blockDim.x * blockIdx.x; + Y0 = threadIdx.y + blockDim.y * blockIdx.y; + // check if we are inside image + if(X0 >= Sz || Y0 >= Sz) return; + // calculate gradient components + int idx = Y0 * Sz + X0; + mirPixSZ *= WEIGHT; + dZdX[idx] = (PAIR(1,0) + DIVSQ2*(PAIR(1,-1)+PAIR(1,1)))/mirPixSZ; + // REMEMBER!!! Axe Y looks up, so we must change the sign + dZdY[idx] = -(PAIR(0,1) + DIVSQ2*(PAIR(1,1)+PAIR(-1,1)))/mirPixSZ; + #undef PAIR + #undef TEX +} + +/** + * Compute matrices of mirror surface deviation + * @param map, mapWH - square matrix of surface deviations (in meters) and its size + * @param mirDia - mirror diameter + * @param mirWH - size of output square matrices (mirWH x mirWH) + * @param mirZ - matrix of mirror Z variations (add to mirror Z) + * @param mirDX, mirDY - partial derivatives dZ/dx & dZ/dy (add to mirror der.) + * @return 0 if fails + */ +EXTERN int CUmkmirDeviat(float *map, size_t mapWH, float mirDia, + mirDeviations * mirDev){ + FNAME(); + ret = 1; + size_t mirWH = mirDev->mirWH; + float *mirDXd = NULL, *mirDYd = NULL, *mirZd = NULL; + size_t mirSp = mirWH*sizeof(float); + size_t mirSz = mirWH*mirSp, dimens = (mirWH+QBLKSZ-1)/QBLKSZ; + size_t pitch; + dim3 blkdim(QBLKSZ, QBLKSZ); + dim3 griddim(dimens, dimens); + // make Z -- simple approximation of + if(!CUbicubic_interpDEV(&mirZd,map,mirWH,mirWH,mapWH,mapWH)){FREERETMACRO;} + CUMOV2HOST(mirDev->mirZ, mirZd, mirSz); + CUFREE(mirZd); + // Z-data would be in pitched 2D array for texture operations + // (to simplify "stretching") + CUALLOCPITCH(mirZd, &pitch, mirSp, mirWH); + CUMOV2DEVPITCH(mirZd, pitch, mirDev->mirZ, mirSp, mirWH); + TEXDATA(mirZd, mirWH, mirWH, pitch); + + CUALLOC(mirDXd, mirSz); + CUALLOC(mirDYd, mirSz); + calcMir_dXdY<<>>(mirWH, mirDia/((float)(mirWH-1)), mirDXd, mirDYd); + cudaThreadSynchronize(); + CUMOV2HOST(mirDev->mirDX, mirDXd, mirSz); + CUMOV2HOST(mirDev->mirDY, mirDYd, mirSz); +free_all: + cudaUnbindTexture(TEXTURE()); + CUFREE(mirDXd); + CUFREE(mirDYd); + CUFREE(mirZd); + return ret; +} +/* + * + * + * <============================================== end of mirror normales block + * + * + */ + +/* + * + * + * Photons trace block ========================================================> + * + * + */ +typedef struct{ + float3 l1; + float3 l2; + float3 l3; +}Matrix; + +// rotation of vec on rotation matrix with lines rM1, rM2, rM3 +__inline__ __device__ float3 rotZA(Matrix *M, float3 vec){ + return make_float3(dot(M->l1,vec), dot(M->l2,vec), dot(M->l3, vec)); +} + +// textures for linear interpolation of mirror deviations +texture TZ; +texture TdX; +texture TdY; +// struct for reducing parameters amount +typedef struct{ + Matrix M; + mirPar P; + BBox B; + bool rot; + Diaphragm D; +}MPB; +/** + * Calculate reflected photon coordinates at plane Z0 + * @param photonX, photonY **INOUT** - + * out: arrays with photon coordinates on image plane + * in: coordinates in diapazone [0,1) + * @param photonSZ - size of prevoius arrays + * @param Z0 - Z-coordinate of image plane + * @param M - mirror rotation matrix or NULL if rotation is absent + * @param mir_WH - mirror derivations matrix size + * @param Parms - mirror parameters + * @param f - light direction vector + * + * @param holes - array with + */ +__global__ void getXY(float *photonX, float *photonY, size_t photonSZ, + float Z0, MPB *mpb, size_t mir_WH, float3 f){ + #define BADPHOTON() do{photonX[IDX] = 1e10f; photonY[IDX] = 1e10f; return;}while(0) + size_t IDX; // pixel number in array + IDX = threadIdx.x + blockDim.x * blockIdx.x; + if(IDX >= photonSZ) return; + Matrix *M = &mpb->M; + mirPar *Parms = &mpb->P; + BBox *box = &mpb->B; + float _2F = Parms->F * 2.f, R = Parms->D / 2.f; + float x,y,z; // coordinates on mirror in meters + x = box->x0 + photonX[IDX] * box->w; // coords of photons + y = box->y0 + photonY[IDX] * box->h; + float r2 = x*x + y*y; + z = r2 / _2F / 2.f; + // we don't mean large inclination so check border by non-inclinated mirror + if(r2 > R*R) BADPHOTON(); + float pixSZ = Parms->D / float(mir_WH-1); // "pixel" size on mirror + // coordinates on deviation matrix, don't forget about y-mirroring! + float xOnMat = (x + R) / pixSZ, yOnMat = (R - y) / pixSZ; + // now add z-deviations, linear interpolation of pre-computed matrix + z += tex2D(TZ, xOnMat, yOnMat); + // point on unrotated mirror + float3 point = make_float3(x,y,z); + // compute normals to unrotated mirror + float3 normal = make_float3( tex2D(TdX, xOnMat, yOnMat) - x/_2F, + tex2D(TdY, xOnMat, yOnMat) - y/_2F, + 1.f + ); + // rotate mirror + if(mpb->rot){ + point = rotZA(M, point); + normal = rotZA(M, normal); + } + normal = normalize(normal); + // calculate reflection direction vector + float3 refl = 2.f*fabs(dot(f, normal))*normal + f; + float K; + if(mpb->D.Nholes){ // there is a diaphragm - test it + Diaphragm *D = &(mpb->D); + int S = D->mask->WH; // size of diaphragm matrix + K = (D->Z - point.z) / refl.z; // scale to convert normal to vector + float xleft = D->box.x0, ybot = D->box.y0; // left bottom angle of dia box + float scalex = D->box.w/(float)S, scaley = D->box.h/(float)S; + x = point.x + K*refl.x; + y = point.y + K*refl.y; + int curX = (int)((x - xleft) / scalex + 0.5f); // coords on dia matrix + int curY = (int)((y - ybot) / scaley + 0.5f); + if(curX < 0 || curX >= S || curY < 0 || curY >= S) BADPHOTON(); + uint16_t mark = D->mask->data[curY*S + curX]; + if(!mark) BADPHOTON(); + do{ + int t = D->holes[mark-1].type; + BBox *b = &D->holes[mark-1].box; + float rx = b->w/2.f, ry = b->h/2.f; + float xc = b->x0 + rx, yc=b->y0 + ry; + float sx = x - xc, sy = y - yc; + switch(t){ + case H_SQUARE: + if(fabs(sx) > rx || fabs(sy) > ry) mark = 0; + break; + case H_ELLIPSE: + if(sx*sx/rx/rx+sy*sy/ry/ry > 1.f) mark = 0; + break; + default: + mark = 0; + } + }while(0); + if(!mark) BADPHOTON(); + } + // OK, test is passed, calculate position on Z0: + K = (Z0 - point.z) / refl.z; + x = point.x + K*refl.x; + y = point.y + K*refl.y; + photonX[IDX] = x; + photonY[IDX] = y; + #undef BADPHOTON +} + +/** + * get coordinates of photons falling onto mirror at point X,Y, reflected + * to plane with z=mirParms->foc + * + * @param xout, yout - arrays of size N_photons with coordinates (allocated outside) + * @param R - use photon coordinates from xout, yout or generate random (in interval [0,1)) + * @param D - deviation & its derivative matrices + * @param mirParms - parameters of mirror + * @param N_photons - number of photons for output (not more than size of input!) + * @param box - bbox where photons are falling + * @return 0 if failed + */ +EXTERN int CUgetPhotonXY(float *xout, float *yout, int R, mirDeviations *D, + mirPar *mirParms, size_t N_photons, BBox *box){ + //FNAME(); + ret = 1; + MPB *mpbdev = NULL, mpb; + aHole *hdev = NULL; + mirMask *mmdev = NULL; + uint16_t *mdatadev = NULL; + float *X = NULL, *Y = NULL; + float SZ = sin(mirParms->objZ); + float z = mirParms->foc; + float *Z_dev = NULL, *dX_dev = NULL, *dY_dev = NULL; + float A = mirParms->Aincl, Z = mirParms->Zincl; + size_t H = D->mirWH, Wp = H*sizeof(float), pitch; + float cA = cos(A), sA = sin(A), cZ = cos(Z), sZ = sin(Z); + size_t sz = N_photons * sizeof(float); + size_t dimens = (N_photons+LBLKSZ-1)/LBLKSZ; + // light direction vector + float3 f = make_float3(-SZ*sin(mirParms->objA), -SZ*cos(mirParms->objA), -cos(mirParms->objZ)); + // rotation matrix by Z than by A: + Matrix M, *Mptr = NULL; + if(A != 0.f || Z != 0.f){ + M.l1 = make_float3(cA, sA*cZ, sA*sZ); + M.l2 = make_float3(-sA, cA*cZ, cA*sZ); + M.l3 = make_float3(0.f, -sZ, cZ); + Mptr = &M; + } + // textures for linear interpolation + CUALLOCPITCH(Z_dev, &pitch, Wp, H); + CUALLOCPITCH(dX_dev, &pitch, Wp, H); + CUALLOCPITCH(dY_dev, &pitch, Wp, H); + CUMOV2DEVPITCH(Z_dev, pitch, D->mirZ, Wp, H); + CUMOV2DEVPITCH(dX_dev, pitch, D->mirDX, Wp, H); + CUMOV2DEVPITCH(dY_dev, pitch, D->mirDY, Wp, H); + #define CUTEX(tex, var) CUTEXTURE(tex, var, H, H, pitch); tex.filterMode = cudaFilterModeLinear + CUTEX(TZ, Z_dev); + CUTEX(TdX, dX_dev); + CUTEX(TdY, dY_dev); + #undef CUTEX + CUALLOC(X, sz); + CUALLOC(Y, sz); + if(R){ + // create __device__ random arrays X & Y + if(!rand_initialized){ + CURAND_CALL(curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT)); + CURAND_CALL(curandSetPseudoRandomGeneratorSeed(gen, throw_random_seed())); + rand_initialized = 1; + } + CURAND_CALL(curandGenerateUniform(gen, X, N_photons)); + cudaThreadSynchronize(); + CURAND_CALL(curandGenerateUniform(gen, Y, N_photons)); + cudaThreadSynchronize(); + }else{ + // move xout, yout to X,Y + CUMOV2DEV(X, xout, sz); + CUMOV2DEV(Y, yout, sz); + } + CUALLOC(mpbdev, sizeof(MPB)); + if(Mptr){ + mpb.rot = true; + memcpy(&mpb.M, Mptr, sizeof(Matrix)); + }else{ + mpb.rot = false; + } + memcpy(&mpb.P, mirParms, sizeof(mirPar)); + memcpy(&mpb.B, box, sizeof(BBox)); + if(!mirParms->dia) mpb.D.Nholes = 0; + // initialize diaphragm + if(mirParms->dia){ // diaphragm is present - allocate memory for it + Diaphragm tmpd; + mirMask tmpM; + memcpy(&tmpd, mirParms->dia, sizeof(Diaphragm)); + memcpy(&tmpM, mirParms->dia->mask, sizeof(mirMask)); + size_t S = sizeof(aHole) * mirParms->dia->Nholes; + CUALLOC(hdev, S); + CUMOV2DEV(hdev, mirParms->dia->holes, S); + tmpd.holes = hdev; + S = mirParms->dia->mask->WH; + S = sizeof(uint16_t) * S * S; + CUALLOC(mdatadev, S); + CUMOV2DEV(mdatadev, mirParms->dia->mask->data, S); + tmpM.data = mdatadev; + CUALLOC(mmdev, sizeof(mirMask)); + CUMOV2DEV(mmdev, &tmpM, sizeof(mirMask)); + tmpd.mask = mmdev; + memcpy(&mpb.D, &tmpd, sizeof(Diaphragm)); + } + CUMOV2DEV(mpbdev, &mpb, sizeof(MPB)); + getXY<<>>(X, Y, N_photons, z, mpbdev, H, f); + cudaThreadSynchronize(); + CUMOV2HOST(xout, X, sz); + CUMOV2HOST(yout, Y, sz); +free_all: + CUFREE(hdev); CUFREE(mmdev); CUFREE(mdatadev); + //CURAND_CALL(curandDestroyGenerator(gen)); + cudaUnbindTexture(TZ); + cudaUnbindTexture(TdX); + cudaUnbindTexture(TdY); + CUFREE(mpbdev); + CUFREE(Z_dev); CUFREE(dX_dev); CUFREE(dY_dev); + CUFREE(X); CUFREE(Y); + return ret; +} +/* + * + * + * <================================================= end of Photons trace block + * + * + */ + + + + +/* +typedef struct{ + float x0; // (x0,y0) - left lower corner + float y0; + float x1; // (x1,y1) - right upper corner + float y1; + size_t pitch; // data pitch in array !!IN NUMBER OF ELEMENTS!! + float dX; // pixel size, dX = (x1-x0) / [image width - 1 ] + float dY; // pixel size, dY = (y1-y0) / [image height - 1] +}devBox; + +__global__ void fillimage(float *image, float *xx, float *yy, size_t sz, devBox *b){ + size_t IDX = threadIdx.x + blockDim.x * blockIdx.x, X, Y; + if(IDX >= sz) return; + float x = xx[IDX], y = yy[IDX]; + if(x < b->x0 || x > b->x1) return; + if(y < b->y0 || y > b->y1) return; + X = (size_t)((x - b->x0) / b->dX); + Y = (size_t)((b->y1 - y) / b->dY); + //atomicAdd(&image[Y*b->pitch + X], 1.f) + image[Y*b->pitch + X] += 1.f; +}*/ +/** + * + * @param phX, phY - photons coordinates + * @param ph_sz - number of photons + * @param image - resulting image (photons **adds** to it) + * @param imW, imH - size of resulting image + * @param imbox - bounding box of resulting image + * @return 0 if fails + */ + /* +EXTERN int CUfillImage(float *phX, float *phY, size_t ph_sz, + float *image, size_t imW, size_t imH, BBox *imbox){ + ret = 1; + size_t linsz = ph_sz * sizeof(float); + size_t pitch, iWp = imW * sizeof(float); + float *xdev = NULL, *ydev = NULL, *imdev = NULL; + devBox *bdev = NULL, devbox; + size_t dimens = (ph_sz+LBLKSZ-1)/LBLKSZ; + + CUALLOC(xdev, linsz); CUALLOC(ydev, linsz); + CUALLOC(bdev, sizeof(devBox)); + CUALLOCPITCH(imdev, &pitch, iWp, imH); +DBG("pitch: %zd", pitch); + CUMOV2DEVPITCH(imdev, pitch, image, iWp, imH); + CUMOV2DEV(xdev, phX, linsz); + CUMOV2DEV(ydev, phY, linsz); + devbox.x0 = imbox->x0; devbox.y0 = imbox->y0; + devbox.x1 = imbox->x0 + imbox->w; + devbox.y1 = imbox->y0 + imbox->h; + devbox.pitch = pitch / sizeof(float); + devbox.dX = imbox->w / (float)(imW - 1); + devbox.dY = imbox->h / (float)(imH - 1); + CUMOV2DEV(bdev, &devbox, sizeof(devBox)); + fillimage<<>>(imdev, xdev, ydev, ph_sz, bdev); + cudaThreadSynchronize(); + CUMOV2HOSTPITCH(image, imdev, pitch, iWp, imH); +free_all: + CUFREE(xdev); CUFREE(ydev); + CUFREE(imdev); CUFREE(bdev); + return ret; +} +*/ diff --git a/src/README b/src/README new file mode 100644 index 0000000..e81d625 --- /dev/null +++ b/src/README @@ -0,0 +1,13 @@ +1. Многопоточные функции поместить в обертку: если есть куда, пытаться сначала вызвать куду, если не сработало или куды нет -- + выполнять на CPU. Все #pragma omp заключить в "скобки": #ifdef OMP ... #endif + +2. Флаг отладки - EBUG + +3. Русифицировать gettext'ом + +4. cmake -DLOCALEDIR=dir изменяет директорию с локалью (снаружи) + + + +ВСЕ ПРОСТРАНСТВЕННЫЕ КООРДИНАТЫ ПРАВОВИНТОВЫЕ +НА ИЗОБРАЖЕНИЯХ ОСЬ Y Н diff --git a/src/cmdlnopts.c b/src/cmdlnopts.c new file mode 100644 index 0000000..6869398 --- /dev/null +++ b/src/cmdlnopts.c @@ -0,0 +1,299 @@ +/* + * 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 +#include +#include */ +#include "cmdlnopts.h" + + +// global variables for parsing +glob_pars G; +mirPar M; +int help; +// don't rewrite existing files: for saveimg.c +extern int rewrite_ifexists; +/* + * variables for mkHartmann.c + */ +extern char *outpfile; // output file name for saving matrixes +extern int printDebug; // print tabulars with data on screen or to file +extern int save_images;// save intermediate results to images + +// suboptions structure for get_suboptions +// in array last MUST BE {0,0,0} +typedef struct{ + float *val; // pointer to result + char *par; // parameter name (CASE-INSENSITIVE!) + bool isdegr; // == TRUE if parameter is an angle in format "[+-][DDd][MMm][SS.S]" +} suboptions; + +// DEFAULTS +// default global parameters +glob_pars const Gdefault = { + 8, // size of initial array of surface deviations + 100, // size of interpolated S0 + 1000, // resulting image size + 10000, // amount of photons falled to one pixel of S1 by one iteration + 0, // add to mask random numbers + 1e-8, // amplitude of added random noice + IT_FITS,// output image type + NULL, // input deviations file name + NULL // mirror +}; +//default mirror parameters +mirPar const Mdefault = { + 6., // diameter + 24.024, // focus + 0., // inclination from Z axe (radians) + 0., // azimuth of inclination (radians) + 0., // azimuth of object (radians) + 0., // zenith of object (radians) + 23.9, // Z-coordinate of light receiver + NULL // diaphragm +}; + +bool get_mir_par(void *arg, int N); +bool get_imtype(void *arg, int N); + +char MirPar[] = N_("set mirror parameters, arg=[diam=num:foc=num:Zincl=ang:Aincl=ang:Ao=ang:Zo=ang:C=num]\n" \ + "\t\t\tALL DEGREES ARE IN FORMAT [+-][DDd][MMm][SS.S] like -10m13.4 !\n" \ + "\t\tdiam - diameter of mirror\n" \ + "\t\tfoc - mirror focus ratio\n" \ + "\t\tZincl - inclination from Z axe\n" \ + "\t\tAincl - azimuth of inclination\n" \ + "\t\tAobj - azimuth of object\n" \ + "\t\tZobj - zenith of object\n" \ + "\t\tccd - Z-coordinate of light receiver"); + +// name has_arg flag val type argptr help +myoption cmdlnopts[] = { + {"help", 0, NULL, 'h', arg_int, APTR(&help), N_("show this help")}, + {"dev-size",1, NULL, 'd', arg_int, APTR(&G.S_dev), N_("size of initial array of surface deviations")}, + {"int-size",1, NULL, 'i', arg_int, APTR(&G.S_interp), N_("size of interpolated array of surface deviations")}, + {"image-size",1,NULL, 'I', arg_int, APTR(&G.S_image), N_("resulting image size")}, + {"N-photons",1,NULL, 'N', arg_int, APTR(&G.N_phot), N_("amount of photons falled to one pixel of matrix by one iteration")}, + {"mir-parameters",1,NULL,'M', arg_function,APTR(&get_mir_par),MirPar}, + {"add-noice",0, &G.randMask,1, arg_none, NULL, N_("add random noice to mirror surface deviations")}, + {"noice-amp",1, NULL, 'a', arg_float, APTR(&G.randAmp), N_("amplitude of random noice (default: 1e-8)")}, + {"dev-file", 1, NULL, 'F', arg_string, APTR(&G.dev_filename),N_("filename for mirror surface deviations (in microns!)")}, + {"force", 0, &rewrite_ifexists,1,arg_none,NULL, N_("rewrite output file if exists")}, + {"log-file",1, NULL, 'l', arg_string, APTR(&outpfile), N_("save matrices to file arg")}, + {"print-matr",0,&printDebug,1, arg_none, NULL, N_("print matrices on screen")}, + {"save-images",0,&save_images,1,arg_none, NULL, N_("save intermediate results to images")}, + {"image-type",1,NULL, 'T', arg_function,APTR(&get_imtype), N_("image type, arg=[jfpt] (Jpeg, Fits, Png, Tiff)")}, + end_option +}; + +/** + * Safely convert data from string to float + * + * @param num (o) - float number read from string + * @param str (i) - input string + * @return TRUE if success + */ +bool myatof(float *num, const char *str){ + float res; + char *endptr; + assert(str); + res = strtof(str, &endptr); + if(endptr == str || *str == '\0' || *endptr != '\0'){ + WARNX(_("Wrong float number format!")); + return FALSE; + } + *num = res; + return TRUE; +} + +/** + * Convert string "[+-][DDd][MMm][SS.S]" into radians + * + * @param ang (o) - angle in radians or exit with help message + * @param str (i) - string with angle + * @return TRUE if OK + */ +bool get_radians(float *ret, char *str){ + float val = 0., ftmp, sign = 1.; + char *ptr; + assert(str); + switch(*str){ // check sign + case '-': + sign = -1.; + case '+': + str++; + } + if((ptr = strchr(str, 'd'))){ // found DDD.DDd + *ptr = 0; if(!myatof(&ftmp, str)) return FALSE; + ftmp = fabs(ftmp); + if(ftmp > 360.){ + WARNX(_("Degrees should be less than 360")); + return FALSE; + } + val += ftmp; + str = ptr + 1; + } + if((ptr = strchr(str, 'm'))){ // found DDD.DDm + *ptr = 0; if(!myatof(&ftmp, str)) return FALSE; + ftmp = fabs(ftmp); + /*if(ftmp >= 60.){ + WARNX(_("Minutes should be less than 60")); + return FALSE; + }*/ + val += ftmp / 60.; + str = ptr + 1; + } + if(strlen(str)){ // there is something more + if(!myatof(&ftmp, str)) return FALSE; + ftmp = fabs(ftmp); + /*if(ftmp >= 60.){ + WARNX(_("Seconds should be less than 60")); + return FALSE; + }*/ + val += ftmp / 3600.; + } + DBG("Angle: %g degr", val*sign); + *ret = D2R(val * sign); // convert degrees to radians + return TRUE; +} + +/** + * 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; +} + +/** + * Parse string of mirror parameters (--mir-diam=...) + * + * @param arg (i) - string with description + * @param N (i) - number of selected option (unused) + * @return TRUE if success + */ +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; +} + +/** + * Get image type from command line parameter: + * F/f for FITS + * P/p for PNG + * J/j for JPEG + * T/t for TIFF + * @param arg (i) - string with parameters + * @return FALSE if fail + */ +bool get_imtype(void *arg, int N _U_){ + assert(arg); + G.it = 0; + do{ + switch(*((char*)arg)){ + case 'F': case 'f': + G.it |= IT_FITS; + break; + case 'P': case 'p': + G.it |= IT_PNG; + break; + case 'J': case 'j': + G.it |= IT_JPEG; + break; + case 'T' : case 't': + G.it |= IT_TIFF; + break; + default: + WARNX(_("Wrong format of image type: %c"), *((char*)arg)); + return FALSE; + } + }while(*((char*)++arg)); + return TRUE; +} diff --git a/src/diaphragm.c b/src/diaphragm.c new file mode 100644 index 0000000..a517d84 --- /dev/null +++ b/src/diaphragm.c @@ -0,0 +1,298 @@ +/* + * diaphragm.c - read diaphragm parameters from 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. + */ + +#include "mkHartmann.h" +#include "wrapper.h" +#include "diaphragm.h" +#include +#include + +/** + * Get double value from object + * @param jobj (i) - object + * @return double value + */ +double get_jdouble(json_object *jobj){ + enum json_type type = json_object_get_type(jobj); + double val; + switch(type){ + case json_type_double: + val = json_object_get_double(jobj); + break; + case json_type_int: + val = json_object_get_int(jobj); + break; + default: + ERRX(_("Wrong value! Get non-number!\n")); + } + return val; +} + +/** + * Fill double array from json array object + * + * @param jobj (i) - JSON array + * @param getLen (o) - array length or NULL + * @return filled array + */ +double *json_get_array(json_object *jobj, int *getLen){ + enum json_type type; + json_object *jarray = jobj; + int arraylen = json_object_array_length(jarray); + double *arr = calloc(arraylen, sizeof(double)); + int L = 0; + int i; + json_object *jvalue; + for (i=0; i< arraylen; i++){ + jvalue = json_object_array_get_idx(jarray, i); + type = json_object_get_type(jvalue); + if(type == json_type_array){ // nested arrays is error + ERRX(_("Invalid file format! Found nested arrays!\n")); + } + else if (type != json_type_object){ + arr[L++] = get_jdouble(jvalue); + } + else{ // non-numerical data? + ERRX(_("Invalid file format! Non-numerical data in array!\n")); + } + } + if(L == 0) FREE(arr); + if(getLen) *getLen = L; + return arr; +} + +/** + * Read bounding box from JSON object + * + * @param B (o) - output BBox (allocated outside) + * @param jobj (i) - JSON object (array with BBox params) + */ +void json_get_bbox(BBox *B, json_object *jobj){ + double *arr = NULL; + int Len; + assert(B); assert(jobj); + json_object *o = json_object_object_get(jobj, "bbox"); + if(!o) return; + if(!(arr = json_get_array(o, &Len)) || Len != 4){ + ERRX(_("\"bbox\" must contain an array of four doubles!\n")); + } + B->x0 = arr[0]; B->y0 = arr[1]; B->w = arr[2]; B->h = arr[3]; + FREE(arr); +} + +aHole globHole; // global parameters for all holes (in beginning of diafragm JSON) +/** + * Fills aHole object by data in JSON + * @param jobj (i) - JSON data for hole + * @param H (o) - output hole object + */ +void get_obj_params(json_object *jobj, aHole *H){ + double *arr = NULL; + int Len; + enum json_type type; + if(!H){ + ERRX( _("Error: NULL instead of aHole structure!\n")); + } + memcpy(H, &globHole, sizeof(aHole)); // initialize hole by global values + json_object *o = json_object_object_get(jobj, "shape"); + if(o){ + const char *ptr = json_object_get_string(o); + if(strcmp(ptr, "square") == 0) H->type = H_SQUARE; + else if(strcmp(ptr, "round") == 0 || strcmp(ptr, "ellipse") == 0 ) H->type = H_ELLIPSE; + else H->type = H_UNDEF; + } + o = json_object_object_get(jobj, "radius"); + if(o){ + type = json_object_get_type(o); + if(type == json_type_int || type == json_type_double){ // circle / square + double R = json_object_get_double(o); + H->box.w = H->box.h = R * 2.; + }else if(type == json_type_array){ // ellipse / rectangle + if(!(arr = json_get_array(o, &Len)) || Len != 2){ + ERRX(_("\"radius\" array must consist of two doubles!\n")); + } + H->box.w = arr[0] * 2.; H->box.h = arr[1] * 2.; + FREE(arr); + }else{ + ERRX(_("\"radius\" must be a number or an array of two doubles!\n")); + } + } + o = json_object_object_get(jobj, "center"); + if(o){ + if(!(arr = json_get_array(o, &Len)) || Len != 2){ + ERRX(_("\"center\" must contain an array of two doubles!\n")); + } + H->box.x0 = arr[0] - H->box.w/2.; + H->box.y0 = arr[1] - H->box.h/2.; + FREE(arr); + }else{ + json_get_bbox(&H->box, jobj); + /* o = json_object_object_get(jobj, "bbox"); + if(o){ + if(!(arr = json_get_array(o, &Len)) || Len != 4){ + ERRX(_("\"bbox\" must contain an array of four doubles!\n")); + } + H->box.x0 = arr[0]; H->box.y0 = arr[1]; H->box.w = arr[2]; H->box.h = arr[3]; + FREE(arr); + }*/ + } +} + +/** + * Fill array of holes from JSON data + * + * @param jobj (i) - JSON array with holes data + * @param getLen (o) - length of array or NULL + * @return holes array + */ +aHole *json_parse_holesarray(json_object *jobj, int *getLen){ + enum json_type type; + json_object *jarray = jobj; + int arraylen = json_object_array_length(jarray), i; + aHole *H = calloc(arraylen, sizeof(aHole)); + json_object *jvalue; + for (i=0; i < arraylen; i++){ + jvalue = json_object_array_get_idx(jarray, i); + type = json_object_get_type(jvalue); + if(type == json_type_object){ + get_obj_params(jvalue, &H[i]); + }else{ + ERRX(_("Invalid holes array format!\n")); + } + } + if(getLen) *getLen = arraylen; + return H; +} + +#ifdef EBUG +char *gettype(aHole *H){ + char *ret; + switch(H->type){ + case H_SQUARE: + ret = "square"; + break; + case H_ELLIPSE: + ret = "ellipse"; + break; + default: + ret = "undefined"; + } + return ret; +} +#endif + +/** + * Try to mmap a file + * + * @param filename (i) - name of file to mmap + * @return pointer with mmap'ed file or die + */ +char *My_mmap(char *filename, size_t *Mlen){ + int fd; + char *ptr; + 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")); + return ptr; +} + +/** + * Read holes array for diaphragm structure from file + * + * file should + * + * @param filename - name of file + * @return readed structure or NULL if failed + */ +aHole *readHoles(char *filename, Diaphragm *dia){ + char *ptr; + enum json_type type; + json_object *o, *jobj; + int i, HolesNum; + aHole *HolesArray; + BBox Dbox = {100.,100.,-200.,-200.}; // bounding box of diaphragm + size_t Mlen; + if(dia) memset(dia, 0, sizeof(Diaphragm)); + ptr = My_mmap(filename, &Mlen); + jobj = json_tokener_parse(ptr); + get_obj_params(jobj, &globHole); // read global parameters + // now try to find diaphragm bounding box & Z-coordinate + json_get_bbox(&Dbox, jobj); + // check for Z-coordinate + if(dia){ + o = json_object_object_get(jobj, "Z"); + if(!o) o = json_object_object_get(jobj, "maskz"); + if(!o) + ERRX(_("JSON file MUST contain floating point field \"Z\" or \"maskz\" with mask's coordinate")); + dia->Z = get_jdouble(o); // read mask Z + } + o = json_object_object_get(jobj, "holes"); + if(!o) + ERRX(_("Corrupted file: no holes found!")); + type = json_object_get_type(o); + if(type == json_type_object){ // single hole + HolesArray = calloc(1, sizeof(aHole)); + assert(HolesArray); + HolesNum = 1; + get_obj_params(o, HolesArray); + }else{ // array of holes + HolesArray = json_parse_holesarray(o, &HolesNum); + } + if(!HolesArray || HolesNum < 1) + ERRX(_("Didn't find any holes in json file!")); + DBG("Readed %d holes", HolesNum); + // check bbox of diafragm (or make it if none) + float minx=100., miny=100., maxx=-100., maxy=-100.; + for(i = 0; i < HolesNum; i++){ + BBox *B = &HolesArray[i].box; + float L=B->x0, R=L+B->w, D=B->y0, U=D+B->h; + if(minx > L) minx = L; + if(maxx < R) maxx = R; + if(miny > D) miny = D; + if(maxy < U) maxy = U; +#ifdef EBUG + green("Hole %4d:", i); + printf(" type =%9s, bbox = {%7.4f, %7.4f, %7.4f, %7.4f }\n", gettype(&HolesArray[i]), + B->x0, B->y0, B->w, B->h); +#endif + } + float wdth = maxx - minx + 0.1, hght = maxy - miny + 0.1; // width & height of bbox + // now correct bbox (or fill it if it wasn't in JSON) + if(Dbox.x0 > minx) Dbox.x0 = minx; + if(Dbox.y0 > miny) Dbox.y0 = miny; + if(Dbox.w < wdth) Dbox.w = wdth; + if(Dbox.h < hght) Dbox.h = hght; + munmap(ptr, Mlen); + if(dia){ + dia->holes = HolesArray; + dia->Nholes = HolesNum; + memcpy(&dia->box, &Dbox, sizeof(BBox)); + } + return HolesArray; +} + diff --git a/src/include/cmdlnopts.h b/src/include/cmdlnopts.h new file mode 100644 index 0000000..d8cbe51 --- /dev/null +++ b/src/include/cmdlnopts.h @@ -0,0 +1,48 @@ +/* + * 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" +#include "mkHartmann.h" +#include "saveimg.h" + +typedef struct{ + int S_dev; // size of initial array of surface deviations + int S_interp; // size of interpolated S0 + int S_image; // resulting image size + int N_phot; // amount of photons falled to one pixel of S1 by one iteration + int randMask; // add to mask random numbers + float randAmp; // amplitude of added random noice + imtype it; // output image type + char *dev_filename;// input deviations file name + mirPar *Mirror; // mirror parameters +} glob_pars; + +// default parameters +extern glob_pars const Gdefault; +extern mirPar const Mdefault; + +glob_pars *parce_args(int argc, char **argv); + +#endif // __CMDLNOPTS_H__ diff --git a/src/include/cutil_math.h b/src/include/cutil_math.h new file mode 100644 index 0000000..c4f829b --- /dev/null +++ b/src/include/cutil_math.h @@ -0,0 +1,772 @@ +/* + * Copyright 1993-2007 NVIDIA Corporation. All rights reserved. + * + * NOTICE TO USER: + * + * This source code is subject to NVIDIA ownership rights under U.S. and + * international Copyright laws. Users and possessors of this source code + * are hereby granted a nonexclusive, royalty-free license to use this code + * in individual and commercial software. + * + * NVIDIA MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THIS SOURCE + * CODE FOR ANY PURPOSE. IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR + * IMPLIED WARRANTY OF ANY KIND. NVIDIA DISCLAIMS ALL WARRANTIES WITH + * REGARD TO THIS SOURCE CODE, INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY, NONINFRINGEMENT, AND FITNESS FOR A PARTICULAR PURPOSE. + * IN NO EVENT SHALL NVIDIA BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL, + * OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS + * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE + * OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE + * OR PERFORMANCE OF THIS SOURCE CODE. + * + * U.S. Government End Users. This source code is a "commercial item" as + * that term is defined at 48 C.F.R. 2.101 (OCT 1995), consisting of + * "commercial computer software" and "commercial computer software + * documentation" as such terms are used in 48 C.F.R. 12.212 (SEPT 1995) + * and is provided to the U.S. Government only as a commercial end item. + * Consistent with 48 C.F.R.12.212 and 48 C.F.R. 227.7202-1 through + * 227.7202-4 (JUNE 1995), all U.S. Government End Users acquire the + * source code with only those rights set forth herein. + * + * Any use of this source code in individual and commercial software must + * include, in the user documentation and internal comments to the code, + * the above Disclaimer and U.S. Government End Users Notice. + */ + +/* + This file implements common mathematical operations on vector types + (float3, float4 etc.) since these are not provided as standard by CUDA. + + The syntax is modelled on the Cg standard library. +*/ + +// This software contains source code provided by NVIDIA Corporation. +// This header contains some bug fixes by Danny Ruijters, changed locations +// are marked with "Danny". +#ifndef CUTIL_MATH_BUGFIXES_H +#define CUTIL_MATH_BUGFIXES_H + +#ifdef CUTIL_MATH_H //Danny +// Houston, we have a problem!!! +// It seems that the buggy cutil_math.h has already been included, and this +// clashes with this header file. +#error cutil_math_bugfixes.h has to be included before cutil_math.h!!! +#endif +#define CUTIL_MATH_H //make sure that the buggy cutil_math.h will not be included + +#include "cuda_runtime.h" + +//////////////////////////////////////////////////////////////////////////////// +typedef unsigned int uint; +typedef unsigned short ushort; + +#ifndef __CUDACC__ +#include + +inline float fminf(float a, float b) +{ + return a < b ? a : b; +} + +inline float fmaxf(float a, float b) +{ + return a < b ? a : b; +} + +inline int max(int a, int b) +{ + return a > b ? a : b; +} + +inline int min(int a, int b) +{ + return a < b ? a : b; +} +#endif + +// float functions +//////////////////////////////////////////////////////////////////////////////// + +// lerp +inline __device__ __host__ float lerp(float a, float b, float t) +{ + return a + t*(b-a); +} + +// clamp +inline __device__ __host__ float clamp(float f, float a, float b) +{ + return fmaxf(a, fminf(f, b)); +} + +// int2 functions +//////////////////////////////////////////////////////////////////////////////// + +// addition +inline __host__ __device__ int2 operator+(int2 a, int2 b) +{ + return make_int2(a.x + b.x, a.y + b.y); +} +inline __host__ __device__ void operator+=(int2 &a, int2 b) +{ + a.x += b.x; a.y += b.y; +} + +// subtract +inline __host__ __device__ int2 operator-(int2 a, int2 b) +{ + return make_int2(a.x - b.x, a.y - b.y); +} +inline __host__ __device__ void operator-=(int2 &a, int2 b) +{ + a.x -= b.x; a.y -= b.y; +} + +// multiply +inline __host__ __device__ int2 operator*(int2 a, int2 b) +{ + return make_int2(a.x * b.x, a.y * b.y); +} +inline __host__ __device__ int2 operator*(int2 a, int s) +{ + return make_int2(a.x * s, a.y * s); +} +inline __host__ __device__ int2 operator*(int s, int2 a) +{ + return make_int2(a.x * s, a.y * s); +} +inline __host__ __device__ void operator*=(int2 &a, int s) +{ + a.x *= s; a.y *= s; +} + +// float2 functions +//////////////////////////////////////////////////////////////////////////////// + +// additional constructors +inline __host__ __device__ float2 make_float2(float s) +{ + return make_float2(s, s); +} +inline __host__ __device__ float2 make_float2(int2 a) +{ + return make_float2(float(a.x), float(a.y)); +} + +// addition +inline __host__ __device__ float2 operator+(float2 a, float2 b) +{ + return make_float2(a.x + b.x, a.y + b.y); +} +inline __host__ __device__ float2 operator+(float2 a, float b) +{ + return make_float2(a.x + b, a.y + b); +} +inline __host__ __device__ float2 operator+(float a, float2 b) +{ + return make_float2(a + b.x, a + b.y); +} +inline __host__ __device__ void operator+=(float2 &a, float2 b) +{ + a.x += b.x; a.y += b.y; +} +inline __host__ __device__ void operator+=(float2 &a, float b) +{ + a.x += b; a.y += b; +} + +// subtract +inline __host__ __device__ float2 operator-(float2 a, float2 b) +{ + return make_float2(a.x - b.x, a.y - b.y); +} +inline __host__ __device__ float2 operator-(float a, float2 b) +{ + return make_float2(a - b.x, a - b.y); +} +inline __host__ __device__ float2 operator-(float2 a, float b) +{ + return make_float2(a.x - b, a.y - b); +} +inline __host__ __device__ void operator-=(float2 &a, float2 b) +{ + a.x -= b.x; a.y -= b.y; +} +inline __host__ __device__ void operator-=(float2 &a, float b) +{ + a.x -= b; a.y -= b; +} + +// multiply +inline __host__ __device__ float2 operator*(float2 a, float2 b) +{ + return make_float2(a.x * b.x, a.y * b.y); +} +inline __host__ __device__ float2 operator*(float2 a, float s) +{ + return make_float2(a.x * s, a.y * s); +} +inline __host__ __device__ float2 operator*(float s, float2 a) +{ + return make_float2(a.x * s, a.y * s); +} +inline __host__ __device__ void operator*=(float2 &a, float s) +{ + a.x *= s; a.y *= s; +} + +// divide +inline __host__ __device__ float2 operator/(float2 a, float2 b) +{ + return make_float2(a.x / b.x, a.y / b.y); +} +inline __host__ __device__ float2 operator/(float2 a, float s) +{ + float inv = 1.0f / s; + return a * inv; +} +inline __host__ __device__ float2 operator/(float s, float2 a) +{ + return make_float2(s / a.x, s / a.y); +} +inline __host__ __device__ void operator/=(float2 &a, float s) +{ + a.x /= s; a.y /= s; +} + +// lerp +inline __device__ __host__ float2 lerp(float2 a, float2 b, float t) +{ + return a + t*(b-a); +} + +// clamp +inline __device__ __host__ float2 clamp(float2 v, float a, float b) +{ + return make_float2(clamp(v.x, a, b), clamp(v.y, a, b)); +} + +inline __device__ __host__ float2 clamp(float2 v, float2 a, float2 b) +{ + return make_float2(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y)); +} + +// dot product +inline __host__ __device__ float dot(float2 a, float2 b) +{ + return a.x * b.x + a.y * b.y; +} + +// length +inline __host__ __device__ float length(float2 v) +{ + return sqrtf(dot(v, v)); +} + +// normalize +inline __host__ __device__ float2 normalize(float2 v) +{ + float invLen = 1.0f / sqrtf(dot(v, v)); + return v * invLen; +} + +// floor +inline __host__ __device__ float2 floor(const float2 v) +{ + return make_float2(floor(v.x), floor(v.y)); +} + +// reflect +inline __host__ __device__ float2 reflect(float2 i, float2 n) +{ + return i - 2.0f * n * dot(n,i); +} + +// float3 functions +//////////////////////////////////////////////////////////////////////////////// + +// additional constructors +inline __host__ __device__ float3 make_float3(float s) +{ + return make_float3(s, s, s); +} +inline __host__ __device__ float3 make_float3(float2 a) +{ + return make_float3(a.x, a.y, 0.0f); +} +inline __host__ __device__ float3 make_float3(float2 a, float s) +{ + return make_float3(a.x, a.y, s); +} +inline __host__ __device__ float3 make_float3(float4 a) +{ + return make_float3(a.x, a.y, a.z); // discards w +} +inline __host__ __device__ float3 make_float3(int3 a) +{ + return make_float3(float(a.x), float(a.y), float(a.z)); +} + +// min +static __inline__ __host__ __device__ float3 fminf(float3 a, float3 b) +{ + return make_float3(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z)); +} + +// max +static __inline__ __host__ __device__ float3 fmaxf(float3 a, float3 b) +{ + return make_float3(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z)); +} + +// addition +inline __host__ __device__ float3 operator+(float3 a, float3 b) +{ + return make_float3(a.x + b.x, a.y + b.y, a.z + b.z); +} +inline __host__ __device__ float3 operator+(float3 a, float b) +{ + return make_float3(a.x + b, a.y + b, a.z + b); +} +inline __host__ __device__ void operator+=(float3 &a, float3 b) +{ + a.x += b.x; a.y += b.y; a.z += b.z; +} + +// subtract +inline __host__ __device__ float3 operator-(float3 a, float3 b) +{ + return make_float3(a.x - b.x, a.y - b.y, a.z - b.z); +} +inline __host__ __device__ float3 operator-(float3 a, float b) +{ + return make_float3(a.x - b, a.y - b, a.z - b); +} +inline __host__ __device__ void operator-=(float3 &a, float3 b) +{ + a.x -= b.x; a.y -= b.y; a.z -= b.z; +} + +// multiply +inline __host__ __device__ float3 operator*(float3 a, float3 b) +{ + return make_float3(a.x * b.x, a.y * b.y, a.z * b.z); +} +inline __host__ __device__ float3 operator*(float3 a, float s) +{ + return make_float3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ float3 operator*(float s, float3 a) +{ + return make_float3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ void operator*=(float3 &a, float s) +{ + a.x *= s; a.y *= s; a.z *= s; +} + +// divide +inline __host__ __device__ float3 operator/(float3 a, float3 b) +{ + return make_float3(a.x / b.x, a.y / b.y, a.z / b.z); +} +inline __host__ __device__ float3 operator/(float3 a, float s) +{ + float inv = 1.0f / s; + return a * inv; +} +inline __host__ __device__ float3 operator/(float s, float3 a) //Danny +{ +// float inv = 1.0f / s; +// return a * inv; + return make_float3(s / a.x, s / a.y, s / a.z); +} +inline __host__ __device__ void operator/=(float3 &a, float s) +{ + float inv = 1.0f / s; + a *= inv; +} + +// lerp +inline __device__ __host__ float3 lerp(float3 a, float3 b, float t) +{ + return a + t*(b-a); +} + +// clamp +inline __device__ __host__ float3 clamp(float3 v, float a, float b) +{ + return make_float3(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b)); +} + +inline __device__ __host__ float3 clamp(float3 v, float3 a, float3 b) +{ + return make_float3(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z)); +} + +// dot product +inline __host__ __device__ float dot(float3 a, float3 b) +{ + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +// cross product +inline __host__ __device__ float3 cross(float3 a, float3 b) +{ + return make_float3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); +} + +// length +inline __host__ __device__ float length(float3 v) +{ + return sqrtf(dot(v, v)); +} + +// normalize +inline __host__ __device__ float3 normalize(float3 v) +{ + float invLen = 1.0f / sqrtf(dot(v, v)); + return v * invLen; +} + +// floor +inline __host__ __device__ float3 floor(const float3 v) +{ + return make_float3(floor(v.x), floor(v.y), floor(v.z)); +} + +// reflect +inline __host__ __device__ float3 reflect(float3 i, float3 n) +{ + return i - 2.0f * n * dot(n,i); +} + +// float4 functions +//////////////////////////////////////////////////////////////////////////////// + +// additional constructors +inline __host__ __device__ float4 make_float4(float s) +{ + return make_float4(s, s, s, s); +} +inline __host__ __device__ float4 make_float4(float3 a) +{ + return make_float4(a.x, a.y, a.z, 0.0f); +} +inline __host__ __device__ float4 make_float4(float3 a, float w) +{ + return make_float4(a.x, a.y, a.z, w); +} +inline __host__ __device__ float4 make_float4(int4 a) +{ + return make_float4(float(a.x), float(a.y), float(a.z), float(a.w)); +} + +// min +static __inline__ __host__ __device__ float4 fminf(float4 a, float4 b) +{ + return make_float4(fminf(a.x,b.x), fminf(a.y,b.y), fminf(a.z,b.z), fminf(a.w,b.w)); +} + +// max +static __inline__ __host__ __device__ float4 fmaxf(float4 a, float4 b) +{ + return make_float4(fmaxf(a.x,b.x), fmaxf(a.y,b.y), fmaxf(a.z,b.z), fmaxf(a.w,b.w)); +} + +// addition +inline __host__ __device__ float4 operator+(float4 a, float4 b) +{ + return make_float4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); +} +inline __host__ __device__ void operator+=(float4 &a, float4 b) +{ + a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w; +} + +// subtract +inline __host__ __device__ float4 operator-(float4 a, float4 b) +{ + return make_float4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); +} +inline __host__ __device__ void operator-=(float4 &a, float4 b) +{ + a.x -= b.x; a.y -= b.y; a.z -= b.z; a.w -= b.w; +} + +// multiply +inline __host__ __device__ float4 operator*(float4 a, float s) +{ + return make_float4(a.x * s, a.y * s, a.z * s, a.w * s); +} +inline __host__ __device__ float4 operator*(float s, float4 a) +{ + return make_float4(a.x * s, a.y * s, a.z * s, a.w * s); +} +inline __host__ __device__ void operator*=(float4 &a, float s) +{ + a.x *= s; a.y *= s; a.z *= s; a.w *= s; +} + +// divide +inline __host__ __device__ float4 operator/(float4 a, float4 b) +{ + return make_float4(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w); +} +inline __host__ __device__ float4 operator/(float4 a, float s) +{ + float inv = 1.0f / s; + return a * inv; +} +inline __host__ __device__ float4 operator/(float s, float4 a) //Danny +{ +// float inv = 1.0f / s; +// return a * inv; + return make_float4(s / a.x, s / a.y, s / a.z, s / a.w); +} +inline __host__ __device__ void operator/=(float4 &a, float s) +{ + float inv = 1.0f / s; + a *= inv; +} + +// lerp +inline __device__ __host__ float4 lerp(float4 a, float4 b, float t) +{ + return a + t*(b-a); +} + +// clamp +inline __device__ __host__ float4 clamp(float4 v, float a, float b) +{ + return make_float4(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b), clamp(v.w, a, b)); +} + +inline __device__ __host__ float4 clamp(float4 v, float4 a, float4 b) +{ + return make_float4(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z), clamp(v.w, a.w, b.w)); +} + +// dot product +inline __host__ __device__ float dot(float4 a, float4 b) +{ + return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; +} + +// length +inline __host__ __device__ float length(float4 r) +{ + return sqrtf(dot(r, r)); +} + +// normalize +inline __host__ __device__ float4 normalize(float4 v) +{ + float invLen = 1.0f / sqrtf(dot(v, v)); + return v * invLen; +} + +// floor +inline __host__ __device__ float4 floor(const float4 v) +{ + return make_float4(floor(v.x), floor(v.y), floor(v.z), floor(v.w)); +} + +// int3 functions +//////////////////////////////////////////////////////////////////////////////// + +// additional constructors +inline __host__ __device__ int3 make_int3(int s) +{ + return make_int3(s, s, s); +} +inline __host__ __device__ int3 make_int3(float3 a) +{ + return make_int3(int(a.x), int(a.y), int(a.z)); +} + +// min +inline __host__ __device__ int3 min(int3 a, int3 b) +{ + return make_int3(min(a.x,b.x), min(a.y,b.y), min(a.z,b.z)); +} + +// max +inline __host__ __device__ int3 max(int3 a, int3 b) +{ + return make_int3(max(a.x,b.x), max(a.y,b.y), max(a.z,b.z)); +} + +// addition +inline __host__ __device__ int3 operator+(int3 a, int3 b) +{ + return make_int3(a.x + b.x, a.y + b.y, a.z + b.z); +} +inline __host__ __device__ void operator+=(int3 &a, int3 b) +{ + a.x += b.x; a.y += b.y; a.z += b.z; +} + +// subtract +inline __host__ __device__ int3 operator-(int3 a, int3 b) +{ + return make_int3(a.x - b.x, a.y - b.y, a.z - b.z); +} + +inline __host__ __device__ void operator-=(int3 &a, int3 b) +{ + a.x -= b.x; a.y -= b.y; a.z -= b.z; +} + +// multiply +inline __host__ __device__ int3 operator*(int3 a, int3 b) +{ + return make_int3(a.x * b.x, a.y * b.y, a.z * b.z); +} +inline __host__ __device__ int3 operator*(int3 a, int s) +{ + return make_int3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ int3 operator*(int s, int3 a) +{ + return make_int3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ void operator*=(int3 &a, int s) +{ + a.x *= s; a.y *= s; a.z *= s; +} + +// divide +inline __host__ __device__ int3 operator/(int3 a, int3 b) +{ + return make_int3(a.x / b.x, a.y / b.y, a.z / b.z); +} +inline __host__ __device__ int3 operator/(int3 a, int s) +{ + return make_int3(a.x / s, a.y / s, a.z / s); +} +inline __host__ __device__ int3 operator/(int s, int3 a) +{ + return make_int3(a.x / s, a.y / s, a.z / s); +} +inline __host__ __device__ void operator/=(int3 &a, int s) +{ + a.x /= s; a.y /= s; a.z /= s; +} + +// clamp +inline __device__ __host__ int clamp(int f, int a, int b) +{ + return max(a, min(f, b)); +} + +inline __device__ __host__ int3 clamp(int3 v, int a, int b) +{ + return make_int3(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b)); +} + +inline __device__ __host__ int3 clamp(int3 v, int3 a, int3 b) +{ + return make_int3(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z)); +} + + +// uint3 functions +//////////////////////////////////////////////////////////////////////////////// + +// additional constructors +inline __host__ __device__ uint3 make_uint3(uint s) +{ + return make_uint3(s, s, s); +} +inline __host__ __device__ uint3 make_uint3(float3 a) +{ + return make_uint3(uint(a.x), uint(a.y), uint(a.z)); +} + +// min +inline __host__ __device__ uint3 min(uint3 a, uint3 b) +{ + return make_uint3(min(a.x,b.x), min(a.y,b.y), min(a.z,b.z)); +} + +// max +inline __host__ __device__ uint3 max(uint3 a, uint3 b) +{ + return make_uint3(max(a.x,b.x), max(a.y,b.y), max(a.z,b.z)); +} + +// addition +inline __host__ __device__ uint3 operator+(uint3 a, uint3 b) +{ + return make_uint3(a.x + b.x, a.y + b.y, a.z + b.z); +} +inline __host__ __device__ void operator+=(uint3 &a, uint3 b) +{ + a.x += b.x; a.y += b.y; a.z += b.z; +} + +// subtract +inline __host__ __device__ uint3 operator-(uint3 a, uint3 b) +{ + return make_uint3(a.x - b.x, a.y - b.y, a.z - b.z); +} + +inline __host__ __device__ void operator-=(uint3 &a, uint3 b) +{ + a.x -= b.x; a.y -= b.y; a.z -= b.z; +} + +// multiply +inline __host__ __device__ uint3 operator*(uint3 a, uint3 b) +{ + return make_uint3(a.x * b.x, a.y * b.y, a.z * b.z); +} +inline __host__ __device__ uint3 operator*(uint3 a, uint s) +{ + return make_uint3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ uint3 operator*(uint s, uint3 a) +{ + return make_uint3(a.x * s, a.y * s, a.z * s); +} +inline __host__ __device__ void operator*=(uint3 &a, uint s) +{ + a.x *= s; a.y *= s; a.z *= s; +} + +// divide +inline __host__ __device__ uint3 operator/(uint3 a, uint3 b) +{ + return make_uint3(a.x / b.x, a.y / b.y, a.z / b.z); +} +inline __host__ __device__ uint3 operator/(uint3 a, uint s) +{ + return make_uint3(a.x / s, a.y / s, a.z / s); +} +inline __host__ __device__ uint3 operator/(uint s, uint3 a) +{ + return make_uint3(a.x / s, a.y / s, a.z / s); +} +inline __host__ __device__ void operator/=(uint3 &a, uint s) +{ + a.x /= s; a.y /= s; a.z /= s; +} + +// clamp +inline __device__ __host__ uint clamp(uint f, uint a, uint b) +{ + return max(a, min(f, b)); +} + +inline __device__ __host__ uint3 clamp(uint3 v, uint a, uint b) +{ + return make_uint3(clamp(v.x, a, b), clamp(v.y, a, b), clamp(v.z, a, b)); +} + +inline __device__ __host__ uint3 clamp(uint3 v, uint3 a, uint3 b) +{ + return make_uint3(clamp(v.x, a.x, b.x), clamp(v.y, a.y, b.y), clamp(v.z, a.z, b.z)); +} + +#endif diff --git a/src/include/diaphragm.h b/src/include/diaphragm.h new file mode 100644 index 0000000..0b45367 --- /dev/null +++ b/src/include/diaphragm.h @@ -0,0 +1,31 @@ +/* + * diaphragm.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 __DIAPHRAGM_H__ +#define __DIAPHRAGM_H__ + +#include "wrapper.h" + +aHole *readHoles(char *filename, Diaphragm *dia); +char *My_mmap(char *filename, size_t *Mlen); + +#endif // __DIAPHRAGM_H__ diff --git a/src/include/mkHartmann.h b/src/include/mkHartmann.h new file mode 100644 index 0000000..dd8dedd --- /dev/null +++ b/src/include/mkHartmann.h @@ -0,0 +1,187 @@ +/* + * mkHartmann.h - main header file with common definitions + * + * 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 __MKHARTMANN_H__ +#define __MKHARTMANN_H__ + +#ifndef _GNU_SOURCE + #define _GNU_SOURCE +#endif +#include +#include // assert +#include // malloc/free etc +#include // printf etc +#include // int types +#include // vprintf +#include // memset etc +#include // file operations +#include // gettext +#include // LONG_MAX +#include // gettimeofday +#include // open +#include // open +#include // open +#include // floor & other math +#include // err +#include +#include // munmap +#include // time, ctime + +#ifndef FLT_EPSILON + #define FLT_EPSILON 1.19209E-07 +#endif + +#ifndef GETTEXT_PACKAGE + #define GETTEXT_PACKAGE "mkHartmann" +#endif +#ifndef PACKAGE_VERSION + #define PACKAGE_VERSION "0.0.1" +#endif +#ifndef LOCALEDIR + #define LOCALEDIR "/usr/share/locale/" +#endif + +#define _(String) gettext(String) +#define gettext_noop(String) String +#define N_(String) gettext_noop(String) + +#define _U_ __attribute__((__unused__)) + +#ifndef THREAD_NUMBER + #define THREAD_NUMBER 2 +#endif + +#ifdef OMP + #ifndef OMP_NUM_THREADS + #define OMP_NUM_THREADS THREAD_NUMBER + #endif + #define Stringify(x) #x + #define OMP_FOR(x) _Pragma(Stringify(omp parallel for x)) +#else + #define OMP_FOR(x) +#endif // OMP + +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) + +// 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 + +#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) + +#ifndef EXTERN // file wasn't included from CUDA.cu + #define EXTERN extern +#endif + +#define RAD 57.2957795130823 +#define D2R(x) ((x) / RAD) +#define R2D(x) ((x) * RAD) + +// STRUCTURES definition + +// bounding box +typedef struct{ + float x0; // left border + float y0; // lower border + float w; // width + float h; // height +} BBox; + +// mirror deviation parameters +typedef struct{ + float *mirZ; // Z + float *mirDX;// dZ/dX + float *mirDY;// dZ/dY + size_t mirWH;// size: mirWH x mirWH +}mirDeviations; + +// type of hole in diaphragm +typedef enum{ + H_SQUARE // square hole + ,H_ELLIPSE // elliptic hole + ,H_UNDEF // error: can't define type +} HoleType; + +// hole in diaphragm +typedef struct{ + BBox box; // bounding box of hole + int type; // type, in case of round hole borders of box are tangents to hole +} aHole; + +// mirror mask for given diaphragm +typedef struct{ + size_t WH; // size of mask: WH x WH + uint16_t *data; // mask data +}mirMask; + +// diaphragm itself +typedef struct{ + BBox box; // bounding box of diaphragm, must be a little larger of its contents + aHole *holes; // array of holes + int Nholes; // size of array + float Z; // z-coordinate of diaphragm + mirMask *mask; +} Diaphragm; + +// parameters of mirror +typedef struct{ + float D; // diameter + float F; // focus + float Zincl; // inclination from Z axe (radians) + float Aincl; // azimuth of inclination (radians) + float objA; // azimuth of object (radians) + float objZ; // zenith of object (radians) + float foc; // Z-coordinate of light receiver + Diaphragm *dia; // diaphragm or NULL +} mirPar; + +// 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); + + +#endif // __MKHARTMANN_H__ + +/* + * + * + * <===========================================================================> + * + * + */ + diff --git a/src/include/parceargs.h b/src/include/parceargs.h new file mode 100644 index 0000000..e5c2479 --- /dev/null +++ b/src/include/parceargs.h @@ -0,0 +1,104 @@ +/* + * 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 + +#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/src/include/saveimg.h b/src/include/saveimg.h new file mode 100644 index 0000000..794c080 --- /dev/null +++ b/src/include/saveimg.h @@ -0,0 +1,42 @@ +/* + * saveimg.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 __SAVEIMG_H__ +#define __SAVEIMG_H__ + +#include +#include "wrapper.h" + +typedef uint8_t imtype; +enum{ + IT_FITS = (imtype)(1), + IT_PNG = (imtype)(1<<1), + IT_JPEG = (imtype)(1<<2), + IT_TIFF = (imtype)(1<<3) +}; + +int writeimg(char *name, imtype t, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data); + +char *createfilename(char* outfile, char* suffix); + +#endif // __SAVEIMG_H__ diff --git a/src/include/wrapper.h b/src/include/wrapper.h new file mode 100644 index 0000000..021f1b3 --- /dev/null +++ b/src/include/wrapper.h @@ -0,0 +1,144 @@ +/* + * wrapper.h - CPU/GPU wrapper + * + * 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 __WRAPPER_H__ +#define __WRAPPER_H__ + +#ifndef EXTERN + #define EXTERN extern +#endif // EXTERN + +#include "mkHartmann.h" + +static const size_t MB = 1024*1024; // convert bytes to MB + +void noCUDA(); +void tryCUDA(); +void forceCUDA(); +int CUsuccess(); +#ifdef CUDA_FOUND + #define CUDAavailable() 1 +#else + #define CUDAavailable() 0 +#endif + +double dtime(); +void getprops(); +EXTERN long throw_random_seed(); +EXTERN int fillImage(float *phX, float *phY, size_t ph_sz, + float *image, size_t imW, size_t imH, BBox *imbox); + +mirMask *makeDmask(Diaphragm *d, size_t minSz, mirPar *M, mirDeviations *D); +void freeDmask(mirMask *m); + + +#define Fn1(A,B) A(x1) +#define Df1(A,B) A(B x1) +#define Fn2(A,B,C) A(x1, x2) +#define Df2(A,B,C) A(B x1, C x2) +#define Fn3(A,B,C,D) A(x1, x2, x3) +#define Df3(A,B,C,D) A(B x1, C x2, D x3) +#define Fn4(A,B,C,D,E) A(x1, x2, x3, x4) +#define Df4(A,B,C,D,E) A(B x1, C x2, D x3, E x4) +#define Fn5(A,B,C,D,E,F) A(x1, x2, x3, x4, x5) +#define Df5(A,B,C,D,E,F) A(B x1, C x2, D x3, E x4, F x5) +#define Fn6(A,B,C,D,E,F,G) A(x1, x2, x3, x4, x5, x6) +#define Df6(A,B,C,D,E,F,G) A(B x1, C x2, D x3, E x4, F x5, G x6) +#define Fn7(A,B,C,D,E,F,G,H) A(x1, x2, x3, x4, x5, x6, x7) +#define Df7(A,B,C,D,E,F,G,H) A(B x1, C x2, D x3, E x4, F x5, G x6, H x7) +#define Fn8(A,B,C,D,E,F,G,H,I) A(x1, x2, x3, x4, x5, x6, x7, x8) +#define Df8(A,B,C,D,E,F,G,H,I) A(B x1, C x2, D x3, E x4, F x5, G x6, H x7, I x8) +#define Fn9(A,B,C,D,E,F,G,H,I,J) A(x1, x2, x3, x4, x5, x6, x7, x8, x9) +#define Df9(A,B,C,D,E,F,G,H,I,J) A(B x1, C x2, D x3, E x4, F x5, G x6, H x7, I x8, J x9) +#define Fn10(A,B,C,D,E,F,G,H,I,J,K) A(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) +#define Df10(A,B,C,D,E,F,G,H,I,J,K) A(B x1, C x2, D x3, E x4, F x5, G x6, H x7, I x8, J x9, K x10) + +#define DEF(N, ...) int Df ## N(__VA_ARGS__) +#define CONCAT(A, B) A ## B +#define FN(N, ...) Fn ## N(__VA_ARGS__) +#define DF(N, ...) Df ## N(__VA_ARGS__) +#define XFUNC(T, X) CONCAT(T, X) +#define FUNC(T, ...) XFUNC(T, FN(__VA_ARGS__)) +#define DFUNC(T,...) EXTERN int XFUNC(T, DF(__VA_ARGS__)) + +#ifdef WRAPPER_C +// even when using cuda in case of fail CUDA init use CPU +static int Only_CPU = +#ifdef CUDA_FOUND + 0 +#else + 1 +#endif +; + +static int CUnoerr, CUforce = 0; + +EXTERN int CUgetprops(); +EXTERN int CUgetMEM(size_t memsz, size_t *Free, size_t *Total); +EXTERN int CUallocaTest(size_t memsz); + +#ifdef CUDA_FOUND +#define SET_F(...) DEF(__VA_ARGS__){ \ + if(!Only_CPU){ CUnoerr = 1; \ + if(FUNC(CU, __VA_ARGS__)) return 1; \ + else CUnoerr = 0; \ + }if(!CUforce && FUNC(CPU, __VA_ARGS__)) return 1; \ + return 0; \ +} +#else +#define SET_F(...) DEF(__VA_ARGS__){ \ + if(FUNC(CPU, __VA_ARGS__)) return 1; \ + return 0; \ +} +#endif // CUDA_FOUND +#else + #define SET_F(...) +#endif // WRAPPER_C + +#ifdef CPU_C // file included from CPU.c + #define BOTH(...) DFUNC(CPU, __VA_ARGS__); + //#pragma message "CPUC" +#elif defined CUDA_CU //file included from CUDA.cu + #define BOTH(...) DFUNC(CU, __VA_ARGS__); +#elif defined WRAPPER_C // wrapper.c needs names of both wariants + #ifndef CUDA_FOUND + #define BOTH(...) DFUNC(CPU, __VA_ARGS__); + #else + #define BOTH(...) DFUNC(CU, __VA_ARGS__); DFUNC(CPU, __VA_ARGS__); + #endif // CUDA_FOUND +#else // file included from something else - just define a function + #define BOTH(...) DFUNC(, __VA_ARGS__); +#endif + +#define DFN(...) BOTH(__VA_ARGS__) SET_F(__VA_ARGS__) + + +DFN(3, fillrandarr, size_t, float *, float) + +DFN(6, bicubic_interp, float *, float *, size_t, size_t, size_t, size_t) + +DFN(4, mkmirDeviat, float *, size_t, float, mirDeviations *) + +DFN(7, getPhotonXY, float *, float *, int, mirDeviations *, mirPar *, size_t, BBox *) + +//DFN(7, fillImage, float *, float *, size_t, float *, size_t, size_t, BBox *) +#endif // __WRAPPER_H__ diff --git a/src/locale/ru/LC_MESSAGES/mkHartmann.mo b/src/locale/ru/LC_MESSAGES/mkHartmann.mo new file mode 100644 index 0000000000000000000000000000000000000000..472b73afe3062ff1a2b5b7be17684936206798df GIT binary patch literal 1118 zcmZ9KPiqrF7{d-qjjNlTFgt1OC-C6e zPvON6;K9@^xV6nD-Ne)i1@$P!77B_d-_0M|I`HsgXP$Z9_ho1G=-$@^mx4NlvQdXo z?@-ZNM;%0Mpe~~RqK>1^?;_+Rau9h8`3QL$SwNmczCs>BzD1rvzDH(IKTuJAY!4xq zkVD85NRB*HumVYWs?2+Kr0(M>{!$ zA|Sb>_H7_>+vYhcE$)D1u7oKrK-aYN5b!A|umTseP^7NQOob$_++f|)?ykus z69vm*FeTW3XL%-}D0q&q5+I5yfwhQAkGcXDlbx5mTLH!o`6N z$ix%(QE`~s^i=MB{MdzCIo9Yq~4>F8E{3@uc}&1 z)pXFW8G7uZ7Sl8|s#!K|Md3Sqx~j)OA24D!JH8Zkr8+LC~{i1Y7Nl_ep6(1d*h*}>^rS95RKNtBej{(_|>tc$M z9!V#+#?|{ZCBYrZ9IQTD#4#j$AuktgYB{$c_mm2k$=lDlHK2+uW7Osp6Kc}Qano}0 z1`Is0q(UHUE-sdrs^vQ1AK3K52DH9Zyali7*P|!2_8n^F#d^bs+EOrI@v6`aYu?vt zppai*{kc!&ZT60_w;$A&yk1Zx2pi2>HJI;)^*ps;p}E>F6}E!;u+f5+U-5!^WvSe= zw7{?WAN=xyzf%6-`|an}M}A}3RX=RE`4MUXw1Snz;Byu0_%PS3cn$DFOkZ9JTII%a nSns^J^>1q%n|}Q6c{TwEUKx literal 0 HcmV?d00001 diff --git a/src/locale/ru/messages.po b/src/locale/ru/messages.po new file mode 100644 index 0000000..3ddd56b --- /dev/null +++ b/src/locale/ru/messages.po @@ -0,0 +1,304 @@ +# 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: 2013-04-02 17:54+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" + +#: mkHartmann.c:168 +#, c-format +msgid "Can't write to %s" +msgstr "" + +#. file not exist but some error occured +#: mkHartmann.c:172 diaphragm.c:216 +#, c-format +msgid "Can't stat %s" +msgstr "" + +#: mkHartmann.c:178 +msgid "No output filename given" +msgstr "" + +#: mkHartmann.c:185 +#, c-format +msgid "Can't open file %s" +msgstr "" + +#: mkHartmann.c:199 +#, c-format +msgid "Can't close file %s" +msgstr "" + +#: mkHartmann.c:229 +msgid "Wrong file: should be matrix of float data separated by spaces" +msgstr "" + +#. update old width counter +#. check it +#: mkHartmann.c:237 +msgid "All rows must contain equal number of columns" +msgstr "" + +#: mkHartmann.c:243 +msgid "Matrix must be square" +msgstr "" + +#: mkHartmann.c:254 +msgid "Input file was modified in runtime!" +msgstr "" + +#: mkHartmann.c:259 +#, c-format +msgid "Error reading data: read %d numbers instaed of %d" +msgstr "" + +#. / "Не могу построить матрицу случайных чисел" +#: mkHartmann.c:316 +msgid "Can't build random matrix" +msgstr "" + +#: mkHartmann.c:323 +msgid "Can't build mirror dZ arrays" +msgstr "" + +#: diaphragm.c:44 +msgid "Wrong value! Get non-number!\n" +msgstr "" + +#. nested arrays is error +#: diaphragm.c:68 +msgid "Invalid file format! Found nested arrays!\n" +msgstr "" + +#. non-numerical data? +#: diaphragm.c:74 +msgid "Invalid file format! Non-numerical data in array!\n" +msgstr "" + +#: diaphragm.c:95 +msgid "\"bbox\" must contain an array of four doubles!\n" +msgstr "" + +#: diaphragm.c:112 +msgid "Error: NULL instead of aHole structure!\n" +msgstr "" + +#: diaphragm.c:130 +msgid "\"radius\" array must consist of two doubles!\n" +msgstr "" + +#: diaphragm.c:135 +msgid "\"radius\" must be a number or an array of two doubles!\n" +msgstr "" + +#: diaphragm.c:141 +msgid "\"center\" must contain an array of two doubles!\n" +msgstr "" + +#: diaphragm.c:178 +msgid "Invalid holes array format!\n" +msgstr "" + +#: diaphragm.c:212 +msgid "No filename given!" +msgstr "" + +#: diaphragm.c:214 +#, c-format +msgid "Can't open %s for reading" +msgstr "" + +#: diaphragm.c:219 +msgid "Mmap error for input" +msgstr "" + +#: diaphragm.c:220 +msgid "Can't close mmap'ed file" +msgstr "" + +#: diaphragm.c:251 +msgid "" +"JSON file MUST contain floating point field \"Z\" or \"maskz\" with mask's " +"coordinate" +msgstr "" + +#: diaphragm.c:256 +msgid "Corrupted file: no holes found!" +msgstr "" + +#: diaphragm.c:267 +msgid "Didn't find any holes in json file!" +msgstr "" + +#. / "Приложение скомпилировано без поддержки CUDA" +#: wrapper.c:39 wrapper.c:48 wrapper.c:58 +msgid "Tool was compiled without CUDA support" +msgstr "" + +#. / "В вычислениях по возможности будет использоваться GPU\n" +#: wrapper.c:74 +msgid "In computations will try to use GPU\n" +msgstr "" + +#. / "Ошибка получения свойств видеоядра" +#: wrapper.c:79 +msgid "Can't get properties" +msgstr "" + +#. / "Ошибка определения доступной памяти" +#: wrapper.c:81 +msgid "Can't determine free memory" +msgstr "" + +#. / "Ошибка выделения памяти" +#: wrapper.c:83 +msgid "Can't allocate memory" +msgstr "" + +#. / "Ошибка в инициализации CUDA!" +#: wrapper.c:99 +msgid "Error in CUDA initialisation!" +msgstr "" + +#. / "ПАМЯТЬ: свободная = " +#: wrapper.c:103 +#, c-format +msgid "MEMORY: free = " +msgstr "" + +#. / " суммарная = " +#: wrapper.c:106 +#, c-format +msgid " total= " +msgstr "" + +#. / "В вычислениях используется только CPU\n" +#: wrapper.c:112 +msgid "Will use only CPU in computations\n" +msgstr "" + +#. / "Тест выделения 100МБ памяти пройден успешно\n" +#: wrapper.c:120 +#, c-format +msgid "Allocation test for 100MB of memory passed\n" +msgstr "" + +#. / "\n\nВсе тесты пройдены успешно" +#: wrapper.c:123 +msgid "" +"\n" +"\n" +"All tests succeed" +msgstr "" + +#. / "Не могу открыть" +#: wrapper.c:154 +msgid "Can't open" +msgstr "" + +#. / "Не могу прочесть" +#: wrapper.c:159 +msgid "Can't read" +msgstr "" + +#: cmdlnopts.c:77 +msgid "" +"set mirror parameters, arg=[diam=num:foc=num:Zincl=ang:Aincl=ang:Ao=ang:" +"Zo=ang:C=num]\n" +"\t\t\tALL DEGREES ARE IN FORMAT [+-][DDd][MMm][SS.S] like -10m13.4 !\n" +"\t\tdiam - diameter of mirror\n" +"\t\tfoc - mirror focus ratio\n" +"\t\tZincl - inclination from Z axe\n" +"\t\tAincl - azimuth of inclination\n" +"\t\tAobj - azimuth of object\n" +"\t\tZobj - zenith of object\n" +"\t\tccd - Z-coordinate of light receiver" +msgstr "" + +#: cmdlnopts.c:89 +msgid "show this help" +msgstr "" + +#: cmdlnopts.c:90 +msgid "size of initial array of surface deviations" +msgstr "" + +#: cmdlnopts.c:91 +msgid "size of interpolated array of surface deviations" +msgstr "" + +#: cmdlnopts.c:92 +msgid "resulting image size" +msgstr "" + +#: cmdlnopts.c:93 +msgid "amount of photons falled to one pixel of matrix by one iteration" +msgstr "" + +#: cmdlnopts.c:95 +msgid "add random noice to mirror surface deviations" +msgstr "" + +#: cmdlnopts.c:96 +msgid "amplitude of random noice (default: 1e-8)" +msgstr "" + +#: cmdlnopts.c:97 +msgid "filename for mirror surface deviations (in microns!)" +msgstr "" + +#: cmdlnopts.c:98 +msgid "rewrite output file if exists" +msgstr "" + +#: cmdlnopts.c:99 +msgid "save matrices to file arg" +msgstr "" + +#: cmdlnopts.c:100 +msgid "print matrices on screen" +msgstr "" + +#: cmdlnopts.c:101 +msgid "save intermediate results to images" +msgstr "" + +#: cmdlnopts.c:102 +msgid "image type, arg=[jfpt] (Jpeg, Fits, Png, Tiff)" +msgstr "" + +#: cmdlnopts.c:119 +msgid "Wrong float number format!" +msgstr "" + +#: cmdlnopts.c:147 +msgid "Degrees should be less than 360" +msgstr "" + +#. wrong format +#: cmdlnopts.c:193 +msgid "Wrong format: no value for keyword" +msgstr "" + +#. nothing found - wrong format +#: cmdlnopts.c:213 +msgid "Bad keyword!" +msgstr "" + +#: cmdlnopts.c:294 +#, c-format +msgid "Wrong format of image type: %c" +msgstr "" diff --git a/src/locale/ru/ru.po b/src/locale/ru/ru.po new file mode 100644 index 0000000..7c669ab --- /dev/null +++ b/src/locale/ru/ru.po @@ -0,0 +1,305 @@ +# 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. +# +msgid "" +msgstr "Project-Id-Version: PACKAGE VERSION\n" + "Report-Msgid-Bugs-To: \n" + "POT-Creation-Date: 2013-04-02 16:58+0400\n" + "PO-Revision-Date: 2013-01-14 18:49+0400\n" + "Last-Translator: Edward V. Emelianov \n" + "Language-Team: LANGUAGE \n" + "Language: Russian\n" + "MIME-Version: 1.0\n" + "Content-Type: text/plain; charset=koi8-r\n" + "Content-Transfer-Encoding: 8bit\n" + +#. / "\n\nВсе тесты пройдены успешно" +#: wrapper.c:123 +#, fuzzy +msgid "\n" + "\n" + "All tests succeed" +msgstr "\n" + "\n" + "Все тесты пройдены успешно\n" + "\n" + +#. / " суммарная = " +#: wrapper.c:106 +#, c-format +msgid " total= " +msgstr " суммарная = " + +#: diaphragm.c:95 +msgid "\"bbox\" must contain an array of four doubles!\n" +msgstr "" + +#: diaphragm.c:141 +msgid "\"center\" must contain an array of two doubles!\n" +msgstr "" + +#: diaphragm.c:130 +msgid "\"radius\" array must consist of two doubles!\n" +msgstr "" + +#: diaphragm.c:135 +msgid "\"radius\" must be a number or an array of two doubles!\n" +msgstr "" + +#. update old width counter +#. check it +#: mkHartmann.c:237 +msgid "All rows must contain equal number of columns" +msgstr "" + +#. / "Тест выделения 100МБ памяти пройден успешно\n" +#: wrapper.c:120 +#, c-format +msgid "Allocation test for 100MB of memory passed\n" +msgstr "Тест выделения 100МБ памяти пройден успешно\n" + +#. nothing found - wrong format +#: cmdlnopts.c:213 +msgid "Bad keyword!" +msgstr "" + +#. / "Ошибка выделения памяти" +#: wrapper.c:83 +msgid "Can't allocate memory" +msgstr "Ошибка выделения памяти" + +#: mkHartmann.c:323 +msgid "Can't build mirror dZ arrays" +msgstr "" + +#. / "Не могу построить матрицу случайных чисел" +#: mkHartmann.c:316 +msgid "Can't build random matrix" +msgstr "" + +#: mkHartmann.c:199 +#, c-format +msgid "Can't close file %s" +msgstr "" + +#: diaphragm.c:220 +msgid "Can't close mmap'ed file" +msgstr "" + +#. / "Ошибка определения доступной памяти" +#: wrapper.c:81 +msgid "Can't determine free memory" +msgstr "Ошибка определения доступной памяти" + +#. / "Ошибка получения свойств видеоядра" +#: wrapper.c:79 +msgid "Can't get properties" +msgstr "Ошибка получения свойств видеоядра" + +#. / "Не могу открыть" +#: wrapper.c:154 +#, fuzzy +msgid "Can't open" +msgstr "Ошибка получения свойств видеоядра" + +#: diaphragm.c:214 +#, c-format +msgid "Can't open %s for reading" +msgstr "" + +#: mkHartmann.c:185 +#, fuzzy, c-format +msgid "Can't open file %s" +msgstr "Ошибка получения свойств видеоядра" + +#. / "Не могу прочесть" +#: wrapper.c:159 +msgid "Can't read" +msgstr "" + +#. file not exist but some error occured +#: mkHartmann.c:172 diaphragm.c:216 +#, c-format +msgid "Can't stat %s" +msgstr "" + +#: mkHartmann.c:168 +#, c-format +msgid "Can't write to %s" +msgstr "" + +#: diaphragm.c:256 +msgid "Corrupted file: no holes found!" +msgstr "" + +#: cmdlnopts.c:147 +msgid "Degrees should be less than 360" +msgstr "" + +#: diaphragm.c:267 +msgid "Didn't find any holes in json file!" +msgstr "" + +#. / "Ошибка в инициализации CUDA!" +#: wrapper.c:99 +msgid "Error in CUDA initialisation!" +msgstr "Ошибка в инициализации CUDA!" + +#: mkHartmann.c:259 +#, c-format +msgid "Error reading data: read %d numbers instaed of %d" +msgstr "" + +#: diaphragm.c:112 +msgid "Error: NULL instead of aHole structure!\n" +msgstr "" + +#. / "В вычислениях по возможности будет использоваться GPU\n" +#: wrapper.c:74 +msgid "In computations will try to use GPU\n" +msgstr "В вычислениях по возможности будет использоваться GPU\n" + +#: mkHartmann.c:254 +msgid "Input file was modified in runtime!" +msgstr "" + +#. nested arrays is error +#: diaphragm.c:68 +msgid "Invalid file format! Found nested arrays!\n" +msgstr "" + +#. non-numerical data? +#: diaphragm.c:74 +msgid "Invalid file format! Non-numerical data in array!\n" +msgstr "" + +#: diaphragm.c:178 +msgid "Invalid holes array format!\n" +msgstr "" + +#: diaphragm.c:251 +msgid "JSON file MUST contain floating point field \"Z\" or \"maskz\" with " + "mask's coordinate" +msgstr "" + +#. / "ПАМЯТЬ: свободная = " +#: wrapper.c:103 +#, c-format +msgid "MEMORY: free = " +msgstr "ПАМЯТЬ: свободная = " + +#: mkHartmann.c:243 +msgid "Matrix must be square" +msgstr "" + +#: diaphragm.c:219 +msgid "Mmap error for input" +msgstr "" + +#: diaphragm.c:212 +msgid "No filename given!" +msgstr "" + +#: mkHartmann.c:178 +msgid "No output filename given" +msgstr "" + +#. / "Приложение скомпилировано без поддержки CUDA" +#: wrapper.c:39 wrapper.c:48 wrapper.c:58 +msgid "Tool was compiled without CUDA support" +msgstr "" + +#. / "В вычислениях используется только CPU\n" +#: wrapper.c:112 +msgid "Will use only CPU in computations\n" +msgstr "В вычислениях используется только CPU\n" + +#: mkHartmann.c:229 +msgid "Wrong file: should be matrix of float data separated by spaces" +msgstr "" + +#: cmdlnopts.c:119 +msgid "Wrong float number format!" +msgstr "" + +#: cmdlnopts.c:294 +#, c-format +msgid "Wrong format of image type: %c" +msgstr "" + +#. wrong format +#: cmdlnopts.c:193 +msgid "Wrong format: no value for keyword" +msgstr "" + +#: diaphragm.c:44 +msgid "Wrong value! Get non-number!\n" +msgstr "" + +#: cmdlnopts.c:95 +msgid "add random noice to mirror surface deviations" +msgstr "" + +#: cmdlnopts.c:93 +msgid "amount of photons falled to one pixel of matrix by one iteration" +msgstr "" + +#: cmdlnopts.c:96 +msgid "amplitude of random noice (default: 1e-8)" +msgstr "" + +#: cmdlnopts.c:97 +msgid "filename for mirror surface deviations (in microns!)" +msgstr "" + +#: cmdlnopts.c:102 +msgid "image type, arg=[jfpt] (Jpeg, Fits, Png, Tiff)" +msgstr "" + +#: cmdlnopts.c:100 +msgid "print matrices on screen" +msgstr "" + +#: cmdlnopts.c:92 +msgid "resulting image size" +msgstr "" + +#: cmdlnopts.c:98 +msgid "rewrite output file if exists" +msgstr "" + +#: cmdlnopts.c:101 +msgid "save intermediate results to images" +msgstr "" + +#: cmdlnopts.c:99 +msgid "save matrices to file arg" +msgstr "" + +#: cmdlnopts.c:77 +msgid "set mirror parameters, arg=[diam=num:foc=num:Zincl=ang:Aincl=ang:" + "Ao=ang:Zo=ang:C=num]\n" + "\t\t\tALL DEGREES ARE IN FORMAT [+-][DDd][MMm][SS.S] like " + "-10m13.4 !\n" + "\t\tdiam - diameter of mirror\n" + "\t\tfoc - mirror focus ratio\n" + "\t\tZincl - inclination from Z axe\n" + "\t\tAincl - azimuth of inclination\n" + "\t\tAobj - azimuth of object\n" + "\t\tZobj - zenith of object\n" + "\t\tccd - Z-coordinate of light receiver" +msgstr "" + +#: cmdlnopts.c:89 +msgid "show this help" +msgstr "" + +#: cmdlnopts.c:90 +msgid "size of initial array of surface deviations" +msgstr "" + +#: cmdlnopts.c:91 +msgid "size of interpolated array of surface deviations" +msgstr "" diff --git a/src/locale/ru/ru.po~ b/src/locale/ru/ru.po~ new file mode 100644 index 0000000..a8febaf --- /dev/null +++ b/src/locale/ru/ru.po~ @@ -0,0 +1,301 @@ +# 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. +# +msgid "" +msgstr "Project-Id-Version: PACKAGE VERSION\n" + "Report-Msgid-Bugs-To: \n" + "POT-Creation-Date: 2013-04-02 16:46+0400\n" + "PO-Revision-Date: 2013-01-14 18:49+0400\n" + "Last-Translator: Edward V. Emelianov \n" + "Language-Team: LANGUAGE \n" + "Language: Russian\n" + "MIME-Version: 1.0\n" + "Content-Type: text/plain; charset=koi8-r\n" + "Content-Transfer-Encoding: 8bit\n" + +#. / "\n\nВсе тесты пройдены успешно" +#: wrapper.c:123 +#, fuzzy +msgid "\n" + "\n" + "All tests succeed" +msgstr "\n" + "\n" + "Все тесты пройдены успешно\n" + "\n" + +#. / " суммарная = " +#: wrapper.c:106 +#, c-format +msgid " total= " +msgstr " суммарная = " + +#: diaphragm.c:95 +msgid "\"bbox\" must contain an array of four doubles!\n" +msgstr "" + +#: diaphragm.c:141 +msgid "\"center\" must contain an array of two doubles!\n" +msgstr "" + +#: diaphragm.c:130 +msgid "\"radius\" array must consist of two doubles!\n" +msgstr "" + +#: diaphragm.c:135 +msgid "\"radius\" must be a number or an array of two doubles!\n" +msgstr "" + +#. update old width counter +#. check it +#: mkHartmann.c:237 +msgid "All rows must contain equal number of columns" +msgstr "" + +#. / "Тест выделения 100МБ памяти пройден успешно\n" +#: wrapper.c:120 +#, c-format +msgid "Allocation test for 100MB of memory passed\n" +msgstr "Тест выделения 100МБ памяти пройден успешно\n" + +#. nothing found - wrong format +#: cmdlnopts.c:213 +msgid "Bad keyword!" +msgstr "" + +#. / "Ошибка выделения памяти" +#: wrapper.c:83 +msgid "Can't allocate memory" +msgstr "Ошибка выделения памяти" + +#. / "Не могу построить матрицу случайных чисел" +#: mkHartmann.c:317 +msgid "Can't build random matrix" +msgstr "" + +#: mkHartmann.c:199 +#, c-format +msgid "Can't close file %s" +msgstr "" + +#: diaphragm.c:220 +msgid "Can't close mmap'ed file" +msgstr "" + +#. / "Ошибка определения доступной памяти" +#: wrapper.c:81 +msgid "Can't determine free memory" +msgstr "Ошибка определения доступной памяти" + +#. / "Ошибка получения свойств видеоядра" +#: wrapper.c:79 +msgid "Can't get properties" +msgstr "Ошибка получения свойств видеоядра" + +#. / "Не могу открыть" +#: wrapper.c:154 +#, fuzzy +msgid "Can't open" +msgstr "Ошибка получения свойств видеоядра" + +#: diaphragm.c:214 +#, c-format +msgid "Can't open %s for reading" +msgstr "" + +#: mkHartmann.c:185 +#, fuzzy, c-format +msgid "Can't open file %s" +msgstr "Ошибка получения свойств видеоядра" + +#. / "Не могу прочесть" +#: wrapper.c:159 +msgid "Can't read" +msgstr "" + +#. file not exist but some error occured +#: mkHartmann.c:172 diaphragm.c:216 +#, c-format +msgid "Can't stat %s" +msgstr "" + +#: mkHartmann.c:168 +#, c-format +msgid "Can't write to %s" +msgstr "" + +#: diaphragm.c:256 +msgid "Corrupted file: no holes found!" +msgstr "" + +#: cmdlnopts.c:147 +msgid "Degrees should be less than 360" +msgstr "" + +#: diaphragm.c:267 +msgid "Didn't find any holes in json file!" +msgstr "" + +#. / "Ошибка в инициализации CUDA!" +#: wrapper.c:99 +msgid "Error in CUDA initialisation!" +msgstr "Ошибка в инициализации CUDA!" + +#: mkHartmann.c:259 +#, c-format +msgid "Error reading data: read %d numbers instaed of %d" +msgstr "" + +#: diaphragm.c:112 +msgid "Error: NULL instead of aHole structure!\n" +msgstr "" + +#. / "В вычислениях по возможности будет использоваться GPU\n" +#: wrapper.c:74 +msgid "In computations will try to use GPU\n" +msgstr "В вычислениях по возможности будет использоваться GPU\n" + +#: mkHartmann.c:254 +msgid "Input file was modified in runtime!" +msgstr "" + +#. nested arrays is error +#: diaphragm.c:68 +msgid "Invalid file format! Found nested arrays!\n" +msgstr "" + +#. non-numerical data? +#: diaphragm.c:74 +msgid "Invalid file format! Non-numerical data in array!\n" +msgstr "" + +#: diaphragm.c:178 +msgid "Invalid holes array format!\n" +msgstr "" + +#: diaphragm.c:251 +msgid "JSON file MUST contain floating point field \"Z\" or \"maskz\" with " + "mask's coordinate" +msgstr "" + +#. / "ПАМЯТЬ: свободная = " +#: wrapper.c:103 +#, c-format +msgid "MEMORY: free = " +msgstr "ПАМЯТЬ: свободная = " + +#: mkHartmann.c:243 +msgid "Matrix must be square" +msgstr "" + +#: diaphragm.c:219 +msgid "Mmap error for input" +msgstr "" + +#: diaphragm.c:212 +msgid "No filename given!" +msgstr "" + +#: mkHartmann.c:178 +msgid "No output filename given" +msgstr "" + +#. / "Приложение скомпилировано без поддержки CUDA" +#: wrapper.c:39 wrapper.c:48 wrapper.c:58 +msgid "Tool was compiled without CUDA support" +msgstr "" + +#. / "В вычислениях используется только CPU\n" +#: wrapper.c:112 +msgid "Will use only CPU in computations\n" +msgstr "В вычислениях используется только CPU\n" + +#: mkHartmann.c:229 +msgid "Wrong file: should be matrix of float data separated by spaces" +msgstr "" + +#: cmdlnopts.c:119 +msgid "Wrong float number format!" +msgstr "" + +#: cmdlnopts.c:294 +#, c-format +msgid "Wrong format of image type: %c" +msgstr "" + +#. wrong format +#: cmdlnopts.c:193 +msgid "Wrong format: no value for keyword" +msgstr "" + +#: diaphragm.c:44 +msgid "Wrong value! Get non-number!\n" +msgstr "" + +#: cmdlnopts.c:95 +msgid "add random noice to mirror surface deviations" +msgstr "" + +#: cmdlnopts.c:93 +msgid "amount of photons falled to one pixel of matrix by one iteration" +msgstr "" + +#: cmdlnopts.c:96 +msgid "amplitude of random noice (default: 1e-8)" +msgstr "" + +#: cmdlnopts.c:97 +msgid "filename for mirror surface deviations (in microns!)" +msgstr "" + +#: cmdlnopts.c:102 +msgid "image type, arg=[jfpt] (Jpeg, Fits, Png, Tiff)" +msgstr "" + +#: cmdlnopts.c:100 +msgid "print matrices on screen" +msgstr "" + +#: cmdlnopts.c:92 +msgid "resulting image size" +msgstr "" + +#: cmdlnopts.c:98 +msgid "rewrite output file if exists" +msgstr "" + +#: cmdlnopts.c:101 +msgid "save intermediate results to images" +msgstr "" + +#: cmdlnopts.c:99 +msgid "save matrices to file arg" +msgstr "" + +#: cmdlnopts.c:77 +msgid "set mirror parameters, arg=[diam=num:foc=num:Zincl=ang:Aincl=ang:" + "Ao=ang:Zo=ang:C=num]\n" + "\t\t\tALL DEGREES ARE IN FORMAT [+-][DDd][MMm][SS.S] like " + "-10m13.4 !\n" + "\t\tdiam - diameter of mirror\n" + "\t\tfoc - mirror focus ratio\n" + "\t\tZincl - inclination from Z axe\n" + "\t\tAincl - azimuth of inclination\n" + "\t\tAobj - azimuth of object\n" + "\t\tZobj - zenith of object\n" + "\t\tccd - Z-coordinate of light receiver" +msgstr "" + +#: cmdlnopts.c:89 +msgid "show this help" +msgstr "" + +#: cmdlnopts.c:90 +msgid "size of initial array of surface deviations" +msgstr "" + +#: cmdlnopts.c:91 +msgid "size of interpolated array of surface deviations" +msgstr "" diff --git a/src/locale/ru/update b/src/locale/ru/update new file mode 100755 index 0000000..968e913 --- /dev/null +++ b/src/locale/ru/update @@ -0,0 +1,2 @@ +#!/bin/sh +msgmerge -Uis ru.po messages.po diff --git a/src/mkHartmann.c b/src/mkHartmann.c new file mode 100644 index 0000000..9e59023 --- /dev/null +++ b/src/mkHartmann.c @@ -0,0 +1,446 @@ +/* + * mkHartmann.c - main file for mkHartmann utilite produsing hartmanogramms + * of BTA telescope + * + * 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 "mkHartmann.h" +#include "cmdlnopts.h" +#include "wrapper.h" +#include "saveimg.h" +#include "diaphragm.h" + +/* + * Coloured messages output + */ +#define RED "\033[1;31;40m" +#define GREEN "\033[1;32;40m" +#define OLDCOLOR "\033[0;0;0m" + + +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; +} + +/* + * 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; +} + +char *outpfile = NULL; // filename for data output in octave text format +int printDebug = 0; // print tab +bool firstRun = TRUE; // first run: create new file +/** + * Print tabular on screen (if outpfile == NULL) or to outpfile + * in octave text format + * Function run only if printDebug == TRUE or outpfile != NULL + * + * @param W, H (i) - matrix width & height + * @param data - data to print + * @param mask - format to print (if on screen) + * @param comment - + * @return + */ +void printTAB(size_t W, size_t H, float *data, char *mask, char *comment){ + size_t x,y; + if(!printDebug && !outpfile) return; // don't print debug info if no flag debug &/| outpfile + if(!outpfile){ // simply print to stdout + if(comment) printf("%s\n", comment); + if(mask) printf("All units *1e-6\n"); + for(y = 0; y < H; y++){ + for(x = 0; x < W; x++){ + float d = data[H*(H-y-1) + x]; + mask ? printf(mask, d*1e6) : printf("%6g ", d); + } + printf("\n"); + } + return; + } + // print to file + struct stat statbuf; + FILE *f = NULL; + #define PR(...) do{if(fprintf(f, __VA_ARGS__) < 0) ERR(_("Can't write to %s"), outpfile);}while(0) + if(firstRun){ + if(stat(outpfile, &statbuf)){ + if(ENOENT != errno) // file not exist but some error occured + ERR(_("Can't stat %s"), outpfile); + // OK, file not exist: use its name + }else{ // file exists, create new name + outpfile = createfilename(outpfile, NULL); // create new file name + } + firstRun = FALSE; + if(!outpfile) ERRX(_("No output filename given")); + f = fopen(outpfile, "w"); // create or truncate + time_t T = time(NULL); + PR("# Created by %s, %s\n", __progname, ctime(&T)); // add header + }else{ // simply open file for adding info + f = fopen(outpfile, "a"); + } + if(!f) ERR(_("Can't open file %s"), outpfile); + // print standard octave header for matrix + if(comment) PR("# name: %s\n", comment); + else PR("# name: tmp_%ju\n", (uintmax_t)time(NULL)); // or create temporary name + PR("# type: matrix\n# rows: %zd\n# columns: %zd\n", H, W); // dimentions + // now put out matrix itself (upside down - for octave/matlab) + for(y = 0; y < H; y++){ + for(x = 0; x < W; x++){ + PR(" %g", *data++); + } + PR("\n"); + } + PR("\n\n"); + #undef PR + if(fclose(f)) ERR(_("Can't close file %s"), outpfile); +} + +/** + * Read array with deviations from file + * if filename is NULL it will just generate zeros (Size x Size) + * ALL DEVIATIONS in file are IN MICROMETERS!!! + * + * @param filename (i) - name of file or NULL for zeros + * @param Size (io) - size of output square array + * @return allocated array + */ +float *read_deviations(char *filename, size_t *Size){ + float *ret = NULL; + int W = 0, W0 = 0, H0 = 0, i; + size_t Mlen; + char *Mem = NULL, *endptr, *ptr; + if(!filename){ + assert(Size); + ret = MALLOC(float, (*Size) * (*Size)); // allocate matrix with given size + assert(ret); + return ret; + } + // there was filename given: try to read data from it + Mem = My_mmap(filename, &Mlen); // from diaphragm.c + ptr = Mem; + do{ + errno = 0; + strtof(ptr, &endptr); + if(errno || (endptr == ptr && *ptr)) + ERRX(_("Wrong file: should be matrix of float data separated by spaces")); + W++; + if(endptr >= Mem + Mlen) break; // eptr out of range - EOF? + if(*endptr == '\n'){ + H0++; + ptr = endptr + 1; + if(!W0) W0 = W; // update old width counter + else if(W != W0) // check it + ERRX(_("All rows must contain equal number of columns")); + W = 0; + }else ptr = endptr; + }while(endptr && endptr < Mem + Mlen); + if(W > 1) H0++; // increase number of rows if there's no trailing '\n' in last line + if(W0 != H0) + ERRX(_("Matrix must be square")); + *Size = W0; + DBG("here"); + ret = MALLOC(float, W0*W0); + DBG("there"); + ptr = Mem; + for(i = 0, H0 = 0; H0 < W0; H0++) + for(W = 0; W < W0; W++, i++){ + DBG("%d ", i); + ret[W0*(W0-H0-1) + W] = strtof(ptr, &endptr) * 1e-6; + if(errno || (endptr == ptr && *ptr) || endptr >= Mem + Mlen) + ERRX(_("Input file was modified in runtime!")); + ptr = endptr; + } + W0 *= W0; + if(i != W0) + ERRX(_("Error reading data: read %d numbers instaed of %d"), W-1, W0); + munmap(Mem, Mlen); + return ret; +} + +int save_images = 0; +/* + * N ph per pix of mask TIME, s + * // 29400 non-zero pts in mask + * 10000 50 + * 50000 73 + * 100000 102 + * 200000 160 + * 500000 303 + * 1000000 556 + * 2000000 1128 + * 5000000 2541 + * 10000000 4826 + */ +int main(int argc, char **argv){ + glob_pars *G = NULL; // default parameters see in cmdlnopts.c + mirPar *M = NULL; // default mirror parameters + int x, y _U_; + // 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); + G = parce_args(argc, argv); + M = G->Mirror; + // Run simple initialisation of CUDA and/or memory test + getprops(); + size_t S0 = G->S_dev, S1 = G->S_interp, Sim = G->S_image, N_phot = G->N_phot; + size_t masksize = S1 * S1; + // type of image to save: fits, png, jpeg or tiff + imtype imt = G->it; + // bounding box of mirror + BBox box; + box.x0 = box.y0 = -M->D/2; box.w = box.h = M->D; + float *idata = read_deviations(G->dev_filename, &S0); + printTAB(S0, S0, idata, "%5.2f ", "input_deviations"); + G->S_dev = S0; // update size + // memory allocation + ALLOC(float, mirZ, masksize); // mirror Z coordinate + ALLOC(float, mirDX, masksize); // projections of normale to mirror + ALLOC(float, mirDY, masksize); + if(G->randMask || G->randAmp != Gdefault.randAmp){ // add random numbers to mask + if(!fillrandarr(S0*S0, idata, G->randAmp)) + /// "Не могу построить матрицу случайных чисел" + ERR(_("Can't build random matrix")); + } + // initialize structure of mirror deviations + mirDeviations mD; + mD.mirZ = mirZ; mD.mirDX = mirDX; mD.mirDY = mirDY; + mD.mirWH = S1; + if(!mkmirDeviat(idata, S0, M->D, &mD)) + ERRX(_("Can't build mirror dZ arrays")); + if(save_images) writeimg("interp_deviations", imt, S1, S1, &box, M, mirZ); + printTAB(S1, S1, mirZ, "%5.3f ", "interpolated_deviations"); + printTAB(S1, S1, mirDX, "%5.2f ", "dev_dZdX"); + printTAB(S1, S1, mirDY, "%5.2f ", "dev_dZdY"); + + /*aHole holes[3] = { + {{-0.35, -0.35, 0.1, 0.1}, H_ELLIPSE}, + {{-0.2, 0.2, 0.7, 0.4}, H_ELLIPSE}, + {{-0.1,-0.45,0.6,0.6}, H_ELLIPSE}}; */ + + // Hartmann mask + Diaphragm dia; //{{-0.5, -0.5, 1., 1.}, NULL, 0, 20., NULL}; + mirMask *mask; + readHoles("holes.json", &dia); + #ifdef EBUG + green("Dia: "); + printf("(%g, %g, %g, %g); %d holes, Z = %g\n", dia.box.x0, dia.box.y0, + dia.box.w, dia.box.h, dia.Nholes, dia.Z); + #endif + if(!(mask = makeDmask(&dia, 128, M, &mD))) ERR("Can't make dmask!"); + M->dia = &dia; + if(save_images){ + int S = (int)dia.mask->WH, S2=S*S; + float *dm = MALLOC(float, S2); + for(x=0; xdata[x]; + writeimg("mirror_mask", imt, S, S, NULL, M, dm); + free(dm); + } + // coordinates of photons + ALLOC(float, xout, N_phot); + ALLOC(float, yout, N_phot); + // resulting image + ALLOC(float, image, Sim*Sim); +/* +//for(int i = 0; i < 100; i++){ + box.x0 = -3.; box.y0 = -3.; box.w = 6.; box.h = 6.; + if(!getPhotonXY(xout, yout, 1, &mD, M, N_phot, &box)) + ERR("Can't build photon map"); + box.x0 = -15e-3; box.y0 = -15e-3; box.w = 30e-3; box.h = 30e-3; + //box.x0 = -5e-3; box.y0 = .8365; box.w = 10e-3; box.h = 10e-3; + if(!fillImage(xout, yout, N_phot, image, Sim, Sim, &box)) + ERR("Can't fill output image"); +//} + writeimg("image", imt, Sim, Sim, &box, M, image); + FREE(xout); FREE(yout); FREE(image); +*/ + // CCD bounding box + BBox CCD = {-15e-3, -15e-3, 30e-3, 30e-3}; + for(x = 0; x < 100; x++){ + if(!getPhotonXY(xout, yout, 1, &mD, M, N_phot, &box)) + ERR("Can't build photon map"); + if(!fillImage(xout, yout, N_phot, image, Sim, Sim, &CCD)) + ERR("Can't fill output image"); + } +/* int S = mask->WH; + double R = M->D / 2., scale = M->D / (double)S; + uint16_t *dptr = mask->data; + box.w = box.h = scale; + // check mask's pixels & throw photons to holes + for(y = 0; y < S; y++){ + for(x = 0; x < S; x++, dptr++){ + if(!*dptr) continue; + DBG("x = %d, Y=%d\n", x,y); + box.x0 = -R + scale*(double)x; + box.y0 = -R + scale*(double)y; + if(!getPhotonXY(xout, yout, 1, &mD, M, N_phot, &box)) + ERR("Can't build photon map"); + if(!fillImage(xout, yout, N_phot, image, Sim, Sim, &CCD)) + ERR("Can't fill output image"); + } + } +*/ + writeimg("image", imt, Sim, Sim, &CCD, M, image); + FREE(xout); FREE(yout); FREE(image); + // if rand() is good, amount of photons on image should be 785398 on every 1000000 + //printTAB(Sim, Sim, image, NULL, "\n\nResulting image:"); +/* for(x = 0; x < N_phot; x++) + if(fabs(xout[x]) < M->D/2. && fabs(yout[x]) < M->D/2.) + printf("photon #%4d:\t\t(%g, %g)\n", x, xout[x]*1e6, yout[x]*1e6);*/ + +/* FILE *F = fopen("TESTs", "w"); + if(!F) ERR("Can't open"); + fprintf(F,"S1\tGPU\t\tCPU\n"); + + + for(S1 = 100; ; S1 += S1*drand48()){ + float *odata = my_alloc(S1*S1, sizeof(float)); + double t0; + fprintf(F,"%zd", S1); + forceCUDA(); + */ +/* int x, y; float *ptr = idata; + printf("Original array:\n"); + for(y = 0; y < 5; y++){ + for(x = 0; x < 5; x++){ + *ptr *= 2.; + *ptr += x; + printf("%4.3f ", (*ptr++)); + } + printf("\n"); + } + t0 = dtime(); + if(!bicubic_interp(odata, idata, S1,S1, S0,S0)) fprintf(F,"\tnan"); + else fprintf(F,"\t%g", dtime() - t0); + + printf("Enlarged array:\n"); + ptr = odata; + for(y = 0; y < 20; y++){ + for(x = 0; x < 20; x++) + printf("%4.3f ", (*ptr++)); + printf("\n"); + } + * + noCUDA(); + t0 = dtime(); + /// "Не могу построить интерполяцию" + if(!bicubic_interp(odata, idata, S1,S1, S0,S0)) ERR(_("Can't do interpolation")); + fprintf(F,"\t%g\n", dtime() - t0); + fflush(F); + free(odata); + }*/ + return 0; +} diff --git a/src/parceargs.c b/src/parceargs.c new file mode 100644 index 0000000..d388ab0 --- /dev/null +++ b/src/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"); + 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/src/saveimg.c b/src/saveimg.c new file mode 100644 index 0000000..885245f --- /dev/null +++ b/src/saveimg.c @@ -0,0 +1,423 @@ +/* + * saveimg.c - functions to save data in png and FITS formats + * + * 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 "mkHartmann.h" +#include "cmdlnopts.h" // for flag "-f", which will tell to rewrite existing file +#include "saveimg.h" +#if defined __PNG && __PNG == TRUE + #include +#endif // PNG +#if defined __JPEG && __JPEG == TRUE + #include +#endif // JPEG +#if defined __TIFF && __TIFF == TRUE + #include +#endif // TIFF +#include + +//char *strdup(const char *s); + +#define TRYFITS(f, ...) \ +do{int status = 0; f(__VA_ARGS__, &status); \ + if (status){ ret = 0; \ + fits_report_error(stderr, status); \ + goto returning;} \ +}while(0) +#define WRITEKEY(...) \ +do{ int status = 0; \ + fits_write_key(fp, __VA_ARGS__, &status);\ + if(status) \ + fits_report_error(stderr, status); \ +}while(0) + +int rewrite_ifexists = 0; // don't rewrite existing files +/** + * Create filename as outfile + number + "." + suffix + * number -- any number from 1 to 9999 + * This function simply returns "outfile.suffix" when "-f" option is set + * + * @param outfile - file name + * @param suffix - file suffix + * @return created filename or NULL + */ +char *createfilename(char* outfile, char* suffix){ + FNAME(); + struct stat filestat; + char buff[256], sfx[32]; + if(suffix) snprintf(sfx, 31, ".%s", suffix); + else sfx[0] = 0; // no suffix + if(rewrite_ifexists){ // there was key "-f": simply return copy of filename + if(snprintf(buff, 255, "%s%s", outfile, sfx) < 1){ + DBG("error"); + return NULL; + } + DBG("(-f present) filename: %s", buff); + return strdup(buff); + } + int num; + if(!outfile) outfile = ""; + for(num = 1; num < 10000; num++){ + if(snprintf(buff, 255, "%s_%04d%s", outfile, num, sfx) < 1){ + DBG("error"); + return NULL; + } + if(stat(buff, &filestat)){ // || !S_ISREG(filestat.st_mode)) // OK, file not exists + DBG("filename: %s", buff); + return strdup(buff); + } + } + DBG("n: %s\n", buff); + WARN("Oops! All numbers are busy or other error!"); + return NULL; +} + +typedef struct{ + float *image; + double min; + double max; + double avr; + double std; +} ImStat; + +static ImStat glob_stat; + +/** + * compute basics image statictics + * @param img - image data + * @param size - image size W*H + * @return + */ +void get_stat(float *img, size_t size){ + FNAME(); + if(glob_stat.image == img) return; + size_t i; + double pv, sum=0., sum2=0., sz=(double)size; + double max = -1., min = 1e15; + for(i = 0; i < size; i++){ + pv = (double) *img++; + sum += pv; + sum2 += (pv * pv); + if(max < pv) max = pv; + if(min > pv) min = pv; + } + glob_stat.avr = sum/sz; + glob_stat.std = sqrt(fabs(sum2/sz - glob_stat.avr*glob_stat.avr)); + glob_stat.max = max; + glob_stat.min = min; + DBG("Image stat: max=%g, min=%g, avr=%g, std=%g", max, min, glob_stat.avr, glob_stat.std); +} + +/** + * Save data to fits file + * @param filename - filename to save to + * @param width, height - image size + * @param imbox - image bounding box + * @data image data + * @return 0 if failed + */ +int writefits(char *filename, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data){ + FNAME(); + long naxes[2] = {width, height}; + char buf[80]; + int ret = 1; + double dX, dY; + if(imbox){ + dX = imbox->w / (double)(width - 1); + dY = imbox->h / (double)(height - 1); + } + time_t savetime = time(NULL); + fitsfile *fp; + get_stat(data, width*height); + assert(filename); + char* newname = MALLOC(char, strlen(filename + 2)); + sprintf(newname, "!%s", filename); // say cfitsio that file could be rewritten + TRYFITS(fits_create_file, &fp, newname); + TRYFITS(fits_create_img, fp, FLOAT_IMG, 2, naxes); + // FILE / Input file original name + WRITEKEY(TSTRING, "FILE", filename, "Input file original name"); + free(newname); + WRITEKEY(TSTRING, "DETECTOR", "Hartmann model", "Detector model"); + if(imbox){ + snprintf(buf, 79, "%.2g x %.2g", dX * 1e6, dY * 1e6); + // PXSIZE / pixel size + WRITEKEY(TSTRING, "PXSIZE", buf, "Pixel size in mkm"); + // XPIXELSZ, YPIXELSZ -- the same + WRITEKEY(TDOUBLE, "XPIXELSZ", &dX, "X pixel size in m"); + WRITEKEY(TDOUBLE, "YPIXELSZ", &dY, "Y pixel size in m"); + // LBCX, LBCY / Coordinates of left bottom corner + WRITEKEY(TFLOAT, "LBCX", &imbox->x0, "X of left bottom corner"); + WRITEKEY(TFLOAT, "LBCY", &imbox->y0, "Y of left bottom corner"); + } + // IMAGETYP / object, flat, dark, bias, scan, eta, neon, push + WRITEKEY(TSTRING, "IMAGETYP", "object", "Image type"); + // DATAMAX, DATAMIN / Max,min pixel value + WRITEKEY(TDOUBLE, "DATAMAX", &glob_stat.max, "Max data value"); + WRITEKEY(TDOUBLE, "DATAMIN", &glob_stat.min, "Min data value"); + // Some Statistics + WRITEKEY(TDOUBLE, "DATAAVR", &glob_stat.avr, "Average data value"); + WRITEKEY(TDOUBLE, "DATASTD", &glob_stat.std, "Standart deviation of data value"); + // DATE / Creation date (YYYY-MM-DDThh:mm:ss, UTC) + strftime(buf, 79, "%Y-%m-%dT%H:%M:%S", gmtime(&savetime)); + WRITEKEY(TSTRING, "DATE", buf, "Creation date (YYYY-MM-DDThh:mm:ss, UTC)"); + // DATE-OBS / DATE OF OBS. + WRITEKEY(TSTRING, "DATE-OBS", buf, "DATE OF OBS. (YYYY-MM-DDThh:mm:ss, local)"); + // OBJECT / Object name + WRITEKEY(TSTRING, "OBJECT", "Modeled object", "Object name"); + // BINNING / Binning + WRITEKEY(TSTRING, "XBIN", "1", "Horizontal binning"); + WRITEKEY(TSTRING, "YBIN", "1", "Vertical binning"); + // PROG-ID / Observation program identifier + WRITEKEY(TSTRING, "PROG-ID", "BTA Hartmann modeling", "Observation program identifier"); + // AUTHOR / Author of the program + WRITEKEY(TSTRING, "AUTHOR", "Edward V. Emelianov", "Author of the program"); + if(mirror){ + WRITEKEY(TFLOAT, "MIRDIAM", &mirror->D, "Mirror diameter"); + WRITEKEY(TFLOAT, "MIRFOC", &mirror->F, "Mirror focus ratio"); + WRITEKEY(TFLOAT, "MIRZINCL", &mirror->Zincl,"Mirror inclination from Z axe"); + WRITEKEY(TFLOAT, "MIRAINCL", &mirror->Aincl,"Azimuth of mirror inclination"); + WRITEKEY(TFLOAT, "A", &mirror->objA, "Object's azimuth"); + WRITEKEY(TFLOAT, "Z", &mirror->objZ, "Object's zenith distance"); + WRITEKEY(TFLOAT, "FOCUS", &mirror->foc, "Z-coordinate of light receiver"); + } + TRYFITS(fits_write_img, fp, TFLOAT, 1, width * height, data); + TRYFITS(fits_close_file, fp); + +returning: + return ret; +} + +static uint8_t *rowptr = NULL; +uint8_t *processRow(float *irow, size_t width, float min, float wd){ + FREE(rowptr); + //float umax = ((float)(UINT16_MAX-1)); + rowptr = MALLOC(uint8_t, width * 3); + OMP_FOR() + for(size_t i = 0; i < width; i++){ + //*ptr = (uint16_t)(umax*(*irow - min)/wd); + double gray = ((double)(irow[i] - min))/((double)wd); + if(gray == 0.) continue; + int G = (int)(gray * 4.); + double x = 4.*gray - (double)G; + uint8_t r = 0, g = 0, b = 0; + uint8_t *ptr = &rowptr[i*3]; + switch(G){ + case 0: + g = (uint8_t)(255. * x + 0.5); + b = 255; + break; + case 1: + g = 255; + b = (uint8_t)(255. * (1. - x) + 0.5); + break; + case 2: + r = (uint8_t)(255. * x + 0.5); + g = 255; + break; + case 3: + r = 255; + g = (uint8_t)(255. * (1. - x) + 0.5); + break; + default: + r = 255; + } + ptr[0] = r; ptr[1] = g; ptr[2] = b; + //ptr[0] = ptr[1] = ptr[2] = gray*255; + } + return rowptr; +} + +int writepng(char *filename, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data){ + FNAME(); + int ret = 1; +#if defined __PNG && __PNG == TRUE + FILE *fp = NULL; + png_structp pngptr = NULL; + png_infop infoptr = NULL; + get_stat(data, width*height); + float min = glob_stat.min, wd = glob_stat.max - min; + float *row; + + if ((fp = fopen(filename, "w")) == NULL){ + perror("Can't open png file"); + ret = 0; + goto done; + } + if ((pngptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, + NULL, NULL, NULL)) == NULL){ + perror("Can't create png structure"); + ret = 0; + goto done; + } + if ((infoptr = png_create_info_struct(pngptr)) == NULL){ + perror("Can't create png info structure"); + ret = 0; + goto done; + } + png_init_io(pngptr, fp); + png_set_compression_level(pngptr, 1); + png_set_IHDR(pngptr, infoptr, width, height, 8, PNG_COLOR_TYPE_RGB,//16, PNG_COLOR_TYPE_GRAY, + PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, + PNG_FILTER_TYPE_DEFAULT); + png_write_info(pngptr, infoptr); + png_set_swap(pngptr); + for(row = &data[width*(height-1)]; height > 0; row -= width, height--) + png_write_row(pngptr, (png_bytep)processRow(row, width, min, wd)); + png_write_end(pngptr, infoptr); +done: + if(fp) fclose(fp); + if(pngptr) png_destroy_write_struct(&pngptr, &infoptr); +#else + WARN("Save to PNG doesn't supported"); +#endif // PNG + return ret; +} + +int writejpg(char *filename, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data){ + FNAME(); + int ret = 1; +#if defined __JPEG && __JPEG == TRUE + get_stat(data, width*height); + float min = glob_stat.min, wd = glob_stat.max - min; + float *row; + FILE* outfile = fopen(filename, "w"); + if(!outfile){ + perror("Can't open jpg file"); + ret = 0; + goto done; + } + struct jpeg_compress_struct cinfo; + struct jpeg_error_mgr jerr; + cinfo.err = jpeg_std_error(&jerr); + jpeg_create_compress(&cinfo); + jpeg_stdio_dest(&cinfo, outfile); + cinfo.image_width = width; + cinfo.image_height = height; + cinfo.input_components = 3; + cinfo.in_color_space = JCS_RGB; + jpeg_set_defaults(&cinfo); + jpeg_set_quality (&cinfo, 99, 1); + jpeg_start_compress(&cinfo, 1); + JSAMPROW row_pointer; + for(row = &data[width*(height-1)]; height > 0; row -= width, height--){ + row_pointer = (JSAMPROW)processRow(row, width, min, wd); + jpeg_write_scanlines(&cinfo, &row_pointer, 1); + } + jpeg_finish_compress(&cinfo); +done: + if(outfile) fclose(outfile); +#else + WARN("Save to JPEG doesn't supported"); +#endif // JPEG + return ret; +} + +int writetiff(char *filename, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data){ + FNAME(); + int ret = 1; +#if defined __TIFF && __TIFF == TRUE + get_stat(data, width*height); + float min = glob_stat.min, wd = glob_stat.max - min; + float *row; + TIFF *image = TIFFOpen(filename, "w"); + if(!image){ + perror("Can't open tiff file"); + ret = 0; + goto done; + } + TIFFSetField(image, TIFFTAG_IMAGEWIDTH, width); + TIFFSetField(image, TIFFTAG_IMAGELENGTH, height); + TIFFSetField(image, TIFFTAG_BITSPERSAMPLE, 8); + TIFFSetField(image, TIFFTAG_SAMPLESPERPIXEL, 3); + TIFFSetField(image, TIFFTAG_ROWSPERSTRIP, 1); + TIFFSetField(image, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT); + TIFFSetField(image, TIFFTAG_COMPRESSION, COMPRESSION_DEFLATE); + TIFFSetField(image, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); + TIFFSetField(image, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); + TIFFSetField(image, TIFFTAG_RESOLUTIONUNIT, RESUNIT_NONE); + tstrip_t strip = 0; + for(row = data; strip < height; row += width, strip++){ + TIFFWriteEncodedStrip(image, strip, processRow(row, width, min, wd), width * 3); + //TIFFWriteScanline + } +done: + if(image) TIFFClose(image); +#else + WARN("Save to TIFF doesn't supported"); +#endif // TIFF + return ret; +} + + +typedef struct{ + imtype t; + char *s; + int (*writefn)(char *, size_t, size_t, BBox *, mirPar *, float *); +}itsuff; + +static itsuff suffixes[] = { + {IT_FITS, "fits", writefits}, + {IT_PNG, "png", writepng}, + {IT_JPEG, "jpeg", writejpg}, + {IT_TIFF, "tiff", writetiff}, + {0, NULL, NULL} +}; +/** + * Save data to image file[s] with format t + * @param name - filename prefix or NULL to save to "outXXXX.format" + * @param t - image[s] type[s] + * @param width, height - image size + * @param imbox - image bounding box (for FITS header) + * @param mirror - mirror parameters (for FITS header) + * @param data image data + * @return number of saved images + */ +int writeimg(char *name, imtype t, size_t width, size_t height, BBox *imbox, + mirPar *mirror, float *data){ + FNAME(); + char *filename = NULL, *suffix; + int ret = 0; + itsuff *suf = suffixes; + while(t && suf->t){ + if(!(t & suf->t)){ + suf++; + continue; + } + t ^= suf->t; + suffix = suf->s; + FREE(filename); + if(name) + filename = createfilename(name, suffix); + else + filename = createfilename("out", suffix); + if(!filename){ + fprintf(stderr, "Create file with name %s and suffix %s failed,\n", name, suffix); + perror("can't make filename"); + continue; + } + if(suf->writefn(filename, width, height, imbox, mirror, data)) ret++; + } + FREE(filename); + return ret; +} diff --git a/src/test/capability.cu b/src/test/capability.cu new file mode 100644 index 0000000..1b284e5 --- /dev/null +++ b/src/test/capability.cu @@ -0,0 +1,25 @@ +#include +#include +#include + +int main(int argc, char **argv){ + cudaDeviceProp dP; + //CUdevice dev; + //CUcontext ctx; + if(cudaSuccess != cudaGetDeviceProperties(&dP, 0)) return 0; + /*if(CUDA_SUCCESS != cuDeviceGet(&dev,0)) return 0; + + // create context for program run: + if(CUDA_SUCCESS != cuCtxCreate(&ctx, 0, dev)) return 0; + printf("\nDevice: %s, totalMem=%zd, memPerBlk=%zd,\n", dP.name, dP.totalGlobalMem, dP.sharedMemPerBlock); + printf("warpSZ=%d, TPB=%d, TBDim=%dx%dx%d\n", dP.warpSize, dP.maxThreadsPerBlock, + dP.maxThreadsDim[0],dP.maxThreadsDim[1],dP.maxThreadsDim[2]); + printf("GridSz=%dx%dx%d, MemovrLap=%d, GPUs=%d\n", dP.maxGridSize[0], + dP.maxGridSize[1],dP.maxGridSize[2], + dP.deviceOverlap, dP.multiProcessorCount); + printf("canMAPhostMEM=%d\n", dP.canMapHostMemory); + printf("compute capability");*/ + printf("-arch=sm_%d%d\n", dP.major, dP.minor); + //cuCtxDetach(ctx); + return 0; +} diff --git a/src/wrapper.c b/src/wrapper.c new file mode 100644 index 0000000..0582b4c --- /dev/null +++ b/src/wrapper.c @@ -0,0 +1,373 @@ +/* + * wrapper.c - CPU / GPU wrapper, try to run function on GPU, if fail - on CPU + * + * 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. + */ + +#define WRAPPER_C +#include "wrapper.h" +#ifdef EBUG + #include "saveimg.h" +#endif + +// Functions to manipulate of forced CPU execution ============================> +void noCUDA(){ // force run on CPU + Only_CPU = 1; + CUforce = 0; +} +void tryCUDA(){ // return to default state +#ifdef CUDA_FOUND + Only_CPU = 0; + CUforce = 0; +#else + /// "Приложение скомпилировано без поддержки CUDA" + WARN(_("Tool was compiled without CUDA support")); +#endif +} +int CUsuccess(){ // test if CU faield +#ifdef CUDA_FOUND + if(Only_CPU) return 0; + else return CUnoerr; +#else + /// "Приложение скомпилировано без поддержки CUDA" + WARN(_("Tool was compiled without CUDA support")); + return 0; +#endif +} +void forceCUDA(){ // not run on CPU even if GPU failed +#ifdef CUDA_FOUND + Only_CPU = 0; + CUforce = 1; +#else + /// "Приложение скомпилировано без поддержки CUDA" + WARN(_("Tool was compiled without CUDA support")); +#endif +} + +// Init function ==============================================================> +/* + * Init CUDA context and/or test memory allocation + * name: getprops + */ +void getprops(){ + FNAME(); + size_t mem = 100 * MB; + /// "Тест на выделение как минимум 100МБ памяти\n" + printf("Make a test for allocation at least 100MB memory\n"); +#ifdef CUDA_FOUND + /// "В вычислениях по возможности будет использоваться GPU\n" + red(_("In computations will try to use GPU\n")); + int status = 0; + char *errors[] = { + "", + /// "Ошибка получения свойств видеоядра" + _("Can't get properties"), + /// "Ошибка определения доступной памяти" + _("Can't determine free memory"), + /// "Ошибка выделения памяти" + _("Can't allocate memory") + }; + size_t theFree, theTotal; + do{ + #define _TEST(F, st) if(!F){status = st;break;} + _TEST(CUgetprops(), 1); + _TEST(CUgetMEM(0, &theFree, &theTotal), 2); + // if there is enough memory + if(theFree / 4 > mem) mem = theFree / 4; + // try to allocate GPU memory + _TEST(CUallocaTest(mem), 3); + #undef _TEST + }while(0); + if(status){ + Only_CPU = 1; + /// "Ошибка в инициализации CUDA!" + WARN(_("Error in CUDA initialisation!")); + WARN(errors[status]); + }else{ + /// "ПАМЯТЬ: свободная = " + printf(_("MEMORY: free = ")); + green("%zdMB,", theFree / MB); + /// " суммарная = " + printf(_(" total= ")); + green("%zdMB\n", theTotal / MB); + } +#endif + if(Only_CPU){ + /// "В вычислениях используется только CPU\n" + green(_("Will use only CPU in computations\n")); + } + // at last step - try to allocate main memory + char *ptr = (char *) malloc(mem); + /// "Ошибка выделения памяти" + if(!ptr || !memset(ptr, 0xaa, mem)) ERR("Error allocating memory"); + free(ptr); + /// "Тест выделения 100МБ памяти пройден успешно\n" + printf(_("Allocation test for 100MB of memory passed\n")); + if(!status){ + /// "\n\nВсе тесты пройдены успешно" + green(_("\n\nAll tests succeed")); + printf("\n\n"); + } +} + +// Functions for pseudo-random number generators initialisation ===============> +/* + * Current time in seconds since UNIX epoch + * name: dtime + * @return time in seconds + */ +double dtime(){ + double t; + struct timeval tv; + gettimeofday(&tv, NULL); + t = tv.tv_sec + ((double)tv.tv_usec)/1e6; + return t; +} +/* + * 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){ + /// "Не могу открыть" + WARN("%s /dev/random!", _("Can't open")); + fail = 1; break; + } + if(sizeof(long) != read(fd, &r_ini, sizeof(long))){ + /// "Не могу прочесть" + WARN("%s /dev/random!", _("Can't read")); + 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); +} + +// Build mask ===================================================> +/** + * Make an array which elements identify corresponding hole's box for given + * box (pixel) on mirror + * It try to make mask size of minSz, but if light from some "pixels" will fall + * onto more than one hole, size would be increased 2 times. Iterations of + * increasing mask size would continue till photons from every "pixel" would + * fall only on one hole + * + * Builded mask attached to input diaphragm as d->mask, before attach a freeDmask() + * is called, so BE CAREFULL! First run of makeDmask() should be with d->mask == NULL! + * + * Mask is an image with cartesian coordinates -> it looks like horizontally mirroring! + * + * @param d - diaphragm for mask making + * @param minSz - minimal mask size + * @param M - mirror parameters + * @param D - mirror deviations + * @return allocated mask (MUST BE FREE by freeDmask) or NULL in case of error + */ +mirMask *makeDmask(Diaphragm *d, size_t minSz, mirPar *M, mirDeviations *D){ + FNAME(); + if(minSz > 4096){ + WARNX("Size of matrix (%d^2) too large!", minSz); + return NULL; + } + if(minSz < 4) minSz = 4; + size_t S2; + int x, y, N; + uint16_t *mdata = NULL, *ddata = NULL; // arrays for diaphragm & mask + mirPar mirror; + memcpy(&mirror, M, sizeof(mirPar)); + mirror.foc = d->Z; // preparing for getPhotonXY + DBG("minSz = %zd", minSz); + S2 = minSz * minSz; + // width & height of "pixel" + double scalex = d->box.w/(double)minSz, scaley = d->box.h/(double)minSz; + + ddata = MALLOC(uint16_t, S2); + // fill diafragm mask to identify which hole to check when "photon" + // falls into that region; fill only BBoxes! + for(N = d->Nholes-1; N > -1; N--){ + BBox *b = &(d->holes[N].box); + // coordinates of hole x0, y0 relative to diaphragm x0,y0 + double X0 = (b->x0 - d->box.x0)/scalex, Y0 = (b->y0 - d->box.y0)/scaley; + // coordinates of upper right corner, add 2. to substitute <= to < in next cycles & to avoid border loss + int x1 = (int)(X0 + b->w/scalex + 2.), y1 = (int)(Y0 + b->h/scaley + 2.); + int x0 = (int)X0, y0 = (int)Y0; + //DBG("scalex: %f scaley: %f, x0=%d, y0=%d, x1=%d, y1=%d\n", scalex, scaley, x0,y0,x1,y1); + uint16_t mark = N + 1; + if(y1 - y0 < 1 || x1 - x0 < 1){ // don't allow scale less than 1pix per hole + DBG("Scale: too little"); + FREE(ddata); + return makeDmask(d, minSz*2, M, D); + } + if(y1 > minSz) y1 = minSz; + if(y0 < 0) y0 = 0; + if(x1 > minSz) x1 = minSz; + if(x0 < 0) x0 = 0; + for(y = y0; y < y1; y++){ + uint16_t *ptr = &ddata[y*minSz + x0]; + for(x = x0; x < x1; x++, ptr++){ + if(*ptr && *ptr != mark){ // zone already occupied, make grid smaller + DBG("Ooops, occupied zone (marker=%d, found %d); try minSz = %zd",mark, *ptr, minSz*2); + FREE(ddata); + return makeDmask(d, minSz*2, M, D); + } + *ptr = mark; + } + } + } + // prepare photons + // they "falls" to corners of inner grid of size (minSz-1)^2 + S2 = (minSz-1)*(minSz-1); + ALLOC(float, Xp, S2); + ALLOC(float, Yp, S2); + // prepare coordinates of "photons" + double mY = 0., dist = 1./(double)(minSz - 2); // distance between points on grid 0..1 + float *xptr = Xp, *yptr = Yp; + for(y = 1; y < minSz; y++, mY += dist){ + double mX = 0.; + for(x = 1; x < minSz; x++, mX += dist){ + *yptr++ = mY; + *xptr++ = mX; + } + } + double mdxy = mirror.D/(double)minSz; // scale on mirror + // box for grid + BBox mirB = {-mirror.D/2.+mdxy, -mirror.D/2.+mdxy, mirror.D-2*mdxy, mirror.D-2*mdxy}; + DBG("mirbox: LD/UR = %g, %g, %g, %g",mirB.x0, mirB.y0, mirB.w+mirB.x0, mirB.h+mirB.y0); + if(!getPhotonXY(Xp, Yp, 0, D, &mirror, S2, &mirB)){ + FREE(Xp); FREE(Yp); FREE(ddata); + return NULL; + } + #ifdef EBUG + float minx=1000.,miny=1000.,maxx=-1000., maxy=-1000.; + for(x=0;x Xp[x])minx=Xp[x]; else if(maxx < Xp[x] && Xp[x] < 1e8) maxx = Xp[x]; + if(miny > Yp[x])miny=Yp[x]; else if(maxy < Yp[x] && Yp[x] < 1e8) maxy = Yp[x];} + DBG("minx: %g, maxx: %g, miny: %g, maxy: %g", minx,maxx,miny,maxy); + #endif + float xleft = d->box.x0, ybot = d->box.y0; + // and now fill mirror mask + int S = minSz - 1; + mdata = MALLOC(uint16_t, minSz*minSz); + ALLOC(uint8_t, histo, d->Nholes); + DBG("xleft: %g, scalex: %g", xleft, scalex); + xptr = Xp; yptr = Yp; + for(y = 0; y < S; y++){ + for(x = 0; x < S; x++, xptr++, yptr++){ + if(*xptr > 1e9 || *yptr > 1e9) continue; // miss to mirror + int curX = (int)((*xptr - xleft) / scalex); + if(curX < 0 || curX >= S) continue; // miss + int curY = (int)((*yptr - ybot) / scaley); + if(curY < 0 ||curY >= S) continue; // miss + uint16_t mark = ddata[curY*minSz + curX]; + if(!mark) continue; // no hole + int pix = y*minSz + x; + int err = 0; + histo[mark-1] = 1; // mark hole as proceeded + #define CHECK(X) if(mdata[X]&&mdata[X]!=mark){err=1;}else{mdata[X]=mark;} + CHECK(pix); // pixel to the left and down + CHECK(pix+1); // right and down + CHECK(pix+minSz); // left & up + CHECK(pix+minSz+1); // right & up + #undef CHECK + if(err){ // zone already occupied, make grid smaller + DBG("Ooops, occupied zone; try minSz = %zd", minSz*2); + FREE(Xp); FREE(Yp); FREE(ddata); FREE(mdata); FREE(histo); + return makeDmask(d, minSz*2, M, D); + } + } + } + FREE(Xp); FREE(Yp); FREE(ddata); + // Now chek whether all holes are present + for(x = 0; x < d->Nholes; x++) + if(!histo[x]){ + DBG("Oooops! Missed a hole!"); + FREE(mdata); + return makeDmask(d, minSz*2, M, D); + } +#ifdef EBUG +int totalpts = 0; S2 = minSz*minSz; uint16_t *ptr = mdata; +for(x=0;xWH = minSz; + Mask->data = mdata; + freeDmask(d->mask); + d->mask = Mask; + return Mask; +} + +void freeDmask(mirMask *m){ + if(!m) return; + FREE(m->data); + FREE(m); +} + +// Build image ================================================================> + +/** + * fill matrix image with photons + * @param phX, phY - photons coordinates + * @param ph_sz - number of photons + * @param image - resulting image (photons **adds** to it) + * @param imW, imH - size of resulting image + * @param imbox - bounding box of resulting image + * @return 0 if fails + */ +int fillImage(float *phX, float *phY, size_t ph_sz, + float *image, size_t imW, size_t imH, BBox *imbox){ + FNAME(); + float x0 = imbox->x0, y0 = imbox->y0, x1 = imbox->x0 + imbox->w, y1 = imbox->y0 + imbox->h; + float dX = imbox->w / (float)(imW - 1), dY = imbox->h / (float)(imH - 1), x=0,y=0; + size_t N; + #ifdef EBUG + float sum = 0., miss = 0., Xc = 0., Yc = 0.; + #endif + for(N = 0; N < ph_sz; N++){ + x = phX[N]; y = phY[N]; + size_t X,Y; + if(x < x0 || x > x1 || y < y0 || y > y1){ + #ifdef EBUG + miss += 1.; + #endif + }else{ + X = (size_t)((x - x0) / dX + 0.5); + Y = (size_t)((y1 - y) / dY + 0.5); + image[Y*imW + X] += 1.f; + #ifdef EBUG + sum += 1.; + Xc += x; Yc += y; + #endif + } + } + DBG("Photons on image: %g, missed: %g; TOTAL: %g\ncenter: Xc=%gmm, Yc=%gmm\nPI=%g", + sum,miss, sum+miss, Xc/sum*1000., Yc/sum*1000., 4.*sum/(sum+miss)); + return 1; +}