Intel SPMD (ISPC) with CMake

####################
# add it to your CMakeLists.txt
#################### ISPC

# given an ispc file, append it and the generated .h file to 'outList'
macro( AddIspcSrc outList outListObj srcLocationsList )
  file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_CFG_INTDIR}/CMakeFiles/ispc)

  foreach(ispc_file ${srcLocationsList})
    message("ISPC source file " ${ispc_file})
    set(srcLocation2 ${ispc_file})
    find_file(srcLocation2 NAMES ${ispc_file} PATHS ${CMAKE_CURRENT_SOURCE_DIR} )
    get_source_file_property(srcLocation2 ${srcLocation2} LOCATION)
    string(FIND "${ispc_file}" "/" lastSlash REVERSE) #just get the file name, strip the path
    if(NOT ${lastSlash} STREQUAL "-1" )
      math(EXPR lastSlash ${lastSlash}+1)
      string(SUBSTRING "${ispc_file}" ${lastSlash} -1 objLocation)
    else()
      set( objLocation ${ispc_file} )
    endif()
    string(REPLACE ".ispc" ".obj" objLocation ${objLocation})
    set(objLocation ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_CFG_INTDIR}/CMakeFiles/ispc/${objLocation}) # put the obj in the intermediate dir
    string(REPLACE ".ispc" "_ispc.h" hdrLocation ${ispc_file})

    message("ISPC header " ${hdrLocation})
    message("ISPC object" ${objLocation})
    ADD_CUSTOM_COMMAND(
       OUTPUT ${objLocation}
       COMMAND ispc ${srcLocation2} -o ${objLocation} -h ${hdrLocation} --arch=x86-64 --target=avx1-i32x8
       IMPLICIT_DEPENDS C ${ispc_file}
       DEPENDS ${ispc_file}
       MAIN_DEPENDENCY  ${ispc_file}
       )
    set_source_files_properties( ${objLocation} PROPERTIES GENERATED TRUE EXTERNAL_OBJECT TRUE )
    set_source_files_properties( ${hdrLocation} PROPERTIES GENERATED TRUE EXTERNAL_OBJECT TRUE )
    list( APPEND ${outList} ${srcLocation} ${hdrLocation} )
    list( APPEND ${outListObj} ${objLocation})
  endforeach(ispc_file)
endmacro( AddIspcSrc )

# collect ispc files
file(GLOB_RECURSE ISPC_FILES  ${PROJECT_SOURCE_DIR}/*.ispc)

# call ispc preprocessor, collect *obj files in outListObj
AddIspcSrc( outList outListObj "${ISPC_FILES}")

add_custom_target(ISPC_OBJECTS ALL DEPENDS ${outListObj})


# your exe file
add_executable(${PROJECT_NAME} ${CPP_FILES})

# make a library of all ispc *obj file
add_library( ispc_objs STATIC EXCLUDE_FROM_ALL ${outListObj})
add_dependencies(ispc_objs ISPC_OBJECTS)
set_target_properties(ispc_objs PROPERTIES LINKER_LANGUAGE C)


# link ispc-ed library with your exe
target_link_libraries(${PROJECT_NAME} ispc_objs)

Unit-testing for CUDA with Google C++ Testing Framework

This code is supposed to be situated in /tests folder as a subdirectory of CMake-based building system

##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
###                     CMakeLists.txt


cmake_minimum_required(VERSION 2.8.7 FATAL_ERROR)
project(tests)

# enable ExternalProject CMake module
include(ExternalProject)

# Set default ExternalProject root directory
set_directory_properties(PROPERTIES EP_PREFIX ${CMAKE_BINARY_DIR}/ThirdParty)

ExternalProject_Add(
    googletest
    URL http://googletest.googlecode.com/files/gtest-1.6.0.zip
    # TIMEOUT 10
    # # Force separate output paths for debug and release builds to allow easy
    # # identification of correct lib in subsequent TARGET_LINK_LIBRARIES commands
    # CMAKE_ARGS -DCMAKE_ARCHIVE_OUTPUT_DIRECTORY_DEBUG:PATH=DebugLibs
    #            -DCMAKE_ARCHIVE_OUTPUT_DIRECTORY_RELEASE:PATH=ReleaseLibs
    #            -Dgtest_force_shared_crt=ON
    # Disable install step
    INSTALL_COMMAND ""
    # Wrap download, configure and build steps in a sсript to log output
    LOG_DOWNLOAD ON
    LOG_CONFIGURE ON
    LOG_BUILD ON)
# Specify include dir
ExternalProject_Get_Property(googletest source_dir)
set(GTEST_INCLUDE_DIR ${source_dir}/include)

# Library
ExternalProject_Get_Property(googletest binary_dir)
set(GTEST_LIBRARY_PATH ${binary_dir}/${CMAKE_FIND_LIBRARY_PREFIXES}gtest.a)
set(GTEST_LIBRARY gtest)

add_library(${GTEST_LIBRARY} UNKNOWN IMPORTED)
set_property(TARGET ${GTEST_LIBRARY}
             PROPERTY IMPORTED_LOCATION ${GTEST_LIBRARY_PATH} )

add_dependencies(${GTEST_LIBRARY} googletest)


# Add test executable target
add_executable(maintest ${PROJECT_SOURCE_DIR}/test.cpp)
# Create dependency of MainTest on googletest
add_dependencies(maintest googletest)

cuda_add_library(cudatests ${GPU_LIBRARY_TYPE}
${PROJECT_SOURCE_DIR}/kernels.cu
OPTIONS ${CUDAOPTIONS}
)

# Specify MainTests link libraries
ExternalProject_Get_Property(googletest binary_dir)

target_link_libraries(maintest ${GTEST_LIBRARY} cudatests)

////////////////////////////////////////////
////////////////////////////////////////////
////////////////////////////////////////////
//                            kernels.cu

#include "kernels.cuh"

__global__
void
square_array_kernel(float* a, unsigned int numElements)
{
  // ### implement me ###
  int i = blockIdx.x * blockDim.x + threadIdx.x;

  if (i < numElements)
    a[i] = a[i] * a[i];
}

void
test()
{
  float*             a_host;                // pointer to array in host memory
  const unsigned int numElements = 100000000;      // number of elements in the array
  size_t             size        = numElements * sizeof(float);
  a_host = new float[size];                 // allocate array on host

  float* a_device;

  cudaMalloc((void**)&a_device, size);
  cudaMemcpy(a_device, a_host, size, cudaMemcpyHostToDevice);

  int block_size = 256;
  int grid_size  = (numElements + block_size - 1) / block_size;

  square_array_kernel << < grid_size, block_size >> > (a_device, numElements);

  cudaMemcpy(a_host, a_device, size, cudaMemcpyDeviceToHost);
  cudaFree(a_device);

  delete[] a_host;
}

////////////////////////////////////////////
////////////////////////////////////////////
////////////////////////////////////////////
//                            kernels.cuh


#ifndef _KERNELS_H_
#define _KERNELS_H_

__global__
void
square_array_kernel(float* a, unsigned int numElements);

void
test();

#endif //  _KERNELS_H_


////////////////////////////////////////////
////////////////////////////////////////////
////////////////////////////////////////////
//                            test.cpp 

#include "gtest/gtest.h"

#include <cuda.h>
#include <cuda_runtime.h>

#include "kernels.cuh"

TEST(A, B) {
  test();

  SUCCEED();
}

int
main(int argc, char** argv)
{
  testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}

Resize images with Python

Recently, I needed to resize some photos in order to upload them to the webpage. There was no point in uploading the full-sized photos, so I downscaled them with a ratio 1/3. The script should be started from the directory with your JPEGs. The output is a sequence of files: resized_0.jpg, resized_1.jpg, resized_2.jpg etc

from os import walk
import Image

def imgResize(im):
    div = 3
    width = im.size[0] / div
    height = im.size[1] / div

    im_result = im.resize((width, height), Image.ANTIALIAS)
    # best down-sizing filter

    return im_result


files = []
for (dirpath, dirnames, filenames) in walk("."):
    files.extend(filenames)
    break

for i, f in enumerate(files):
    if "JPG" in f and not "resized" in f:
        im = Image.open(f)
        im = imgResize(im)
        im.save("resized_"+str(i)+".jpg")

Многомерная сплайн-интерполяция. Построение карт

Моя старая магистерская диссертация на тему многомерной интерполяции данных в геофизике.

https://docs.google.com/file/d/0B32kbaZOrKltT0NJSWphUHJySnM/edit

How to make anaglyph (stereo, 3D) image with standard camera

If you do not have special equipment you still can create pretty good anaglyph images with standard camera.

What do we need for it?

1) Make two shots with little parallax. Try not to make large offset. For example if you make photo of a tree then offset must be about 6-10 сm.

2) You need special anaglyph glasses red-blue glasses.



3) Make following steps to get final image:

- Apply screen blenging operation (wiki) to the left and right images. Left must be blended with cyan color (0x00FFFF) and right with red color (0xFF0000).



- Having two changed images for left and right eyes blend them with multiply operation (wiki).



Make sure you adjusted images in such way, that objects on the photos matches with each other.

- Finally you get result. Put on glasses and enjoy!

Трехмерная модель скважины

Эту программу я начал писать еще два года назад как курсовую работу, она так и осталась недоделанная.

Недавно потребовалось довести проект до ума, в результате чего более половины исходников было переписано и добавлены новые возможности.

Программа достает из базы данных данные инклинометрии и любое количество кривых каверномера (профилемера).

Для расчета изменения направления ствола скважины использованы преобразования координат при помощи матриц поворота.

Программа позволяет наглядно представить траекторию и форму ствола скважины, оценить образовавшиеся каверны, провести сопоставление отклонений от номинального диаметра с глубиной залегания нефтеносных пластов.

В результате получаем такие картинки:

Построение изолиний поверхности

Решение задачи не претендует на абюсолютную правильность реализации и не явлеяется самым лучшим и оптимальным. Предложенный алгоритм – попытка решения поставленной проблемы.

Задача.

Поверхность – двумерный массив размером nx на ny точек, в каждой из которых содержится Z отметка высоты. Пусть, поверхность хранится в массиве вида:

std::vector<double> surf;

К элементу с номером i, j можно обратиться следующим образом: surf[j*nx + i];

Задача состоит в трассировке линий, принадлежащих этой поверхности, состоящих из точек с одинаковыми Z-отметками.

Элементарные элементы из которых состоит поверхность – ячейки – прямоугольники, образованные четыремя соседними вершинами. Эти ячейки не подходят для трассировки, т.к. через каждый прямоугольник можно провести две изолинии, притом, несколькими способами.

Треугольник является такой фигурой, через которую можно провести только одну изолинию. Это означает что если плоскость, описываемая уравнением Z = const пересекает какую-либо сторону треугольника, она гарантировано пересечет одну из двух оставшихся сторон. Это легко доказывается.

Поэтому для трассировки изолиний необходимо разбить существующие прямоугольники на треугольники следующим образом:



Подготовка.

Трассировка изолиний выполняется с заданным шагом dz в интервале изменения высот поверхности. Далее приведено описание построения изолинии для одной величины Z.

Первым этапом является проверка всех существующих граней на пересечение с заданной величиной Z.

Каждую грань легко задать одним узлом и направлением. Всего существует три направления: горизонтальное, вертикальное, диагональное. Поэтому, для хранения отмеченных граней понадобится массив такого же размера, как и масив поверхности nx на ny. Каждая ячейка массива будет хранить unsigned char число, являющееся битовым полем. В битовом поле закодированны три направления: 001 – горизонтальное, 010 – вертикальное, 100 – диагональное. Если ячейка массива содержит 0 – значит из этого узла не выходит ни одного ребра, пересекающегося с заданной величиной Z.

std::vector<unsigned char> edges; // bit





// так можно отметить вертикальную грань


if ((z1 < curz) && (z >= curz))


edges[y * nx + x] |= 2;

Итак, в цикле перебираем все точки поверхности и заносим информацию об отмеченных ребрах с специальный массив.



Трассировка.

Для того что бы трассировать линии удобно применить такой контейнер как очередь (FILO – First In Last Out). На каждом шаге цикла из очереди будет доставаться одна грань, рассчитываться координаты точки пересечения этой грани с плоскостью Z = Z0, анализироваться треугольник, которому принадлежит грань и добавляться новые грани в конец очереди.

Элементы хранимые в очереди – грани треугольников – отрезки соединяющие две соседние точки в массиве поверхности. Грань может быть представлена такой структурой данных:

struct Edge


{


int knot1;


int knot2;


};

Каждый узел поверхности однозначно определяется целым числом, значение которого можно вычислить как j*nx+i, где i и j – координаты узла в массиве.
Так как каждая грань принадлежит двум треугольникам, удобно задавать ее координатами начала и конца. Это позволяет определить направление вектора-грани, что дает информацию о треугольнике, к которому принадлежит грань. Будем считать что все треугольники обходятся против часовой стрелки, и что бы узнать третью координату треугольника, достаточно обойти его против часовой стрелки по направлению вектора-грани.

Первый этап

Для начала надо найти хотя бы одну отмеченную грань. Последовательно перебираем массив, где хранятся метки граней и ищем значение не равное 0. Это означает, что из данный узел является началом хотя бы одной грани (горизонтальной, вертикальной, диагональной), пересекающей поверхность Z = Z0;

После того как узел найден, проверяем какая грань, выходящая из него, отмечена и добавляем эту грань в очередь. Причем добавляется прямое и инвертированное направление грани: (knot1, knot2) и (knot2, knot1). Это позволяет начать трассировку линий в двух направлениях.

Втоой этап.

После того, как в очеред внесены первые две грани, необходимо запустить цикл по очереди.

а) Достаем грань из очереди и находим третью вершину, вместе с которой она образует треугольник. Третью вершину легко найти зная координаты вершин грани (knot1, knot2).
б) Проверяем две оставшиеся грани реугольника, если какая-либо грань пересекает плоскоссть Z = Z0 (проверяем массив с отмечеными гранями), добавляем ее в очередь.
в) Считаем точные координаты точки пересечения исходной грани треугольника с горизонтальной плоскостью. Это просто сделать, зная значения высот узлов грани Z1, Z2. и координаты точек (x1, y1) (x2, y2)
г) Добавляем полученную координату к текущей изолинии. Убираем метку этой грани из массива отмеченных граней и удяляем грань из очереди.

Повторяем цикл, пока очередь не исчерпается.

Стоит заметить, что самая первая грань в очереди добавляется два раза, в прямом и инвертированном направлении. Это необходимо для того, что бы начать построение изолинии сразу в двух направлениях. Если изолиния замкнутая, к таком способу можно не прибегать, но если изолиния разомкнутая и заканчивается на границах, требуется отследить ее в обе стороны, пока изолиния не дойдет до границы.. Для того что бы знать, в начало или в конец изолинии добавлять новую точку, можно хранить вместе с каждой отмеченной гранью переменную int direction, которая может принимать значения +1 или -1.

struct Edge


{





int knot1;






int knot2;










int direction;



};

Я использовал слово “грань”, которое в общем случае обозначает поверхность объемного многогранница для удобства, т.к. словосочетание “сторона треугольника” кажется мне не ёмким