dahtah / imager

R package for image processing
GNU Lesser General Public License v3.0
187 stars 43 forks source link

Convert to polar coordinates and back? #137

Open mconsidine opened 3 years ago

mconsidine commented 3 years ago

Hi, I'm dealing with an image of the sun and would like to apply a function to a part of the image around the perimeter of the solar disk. It's occurring to me that life might be easier if I converted the image to polar coordinates, applied whatever transformation I want to and then convert it back.

Before I end up reinventing the wheel, does anyone know if there is some function(s) buried in imager that can help me accomplish this?

EDIT: Someone asked a similar question related to python here https://stackoverflow.com/questions/51675940/converting-an-image-from-cartesian-to-polar-limb-darkening

TIA, mconsidine

mconsidine commented 3 years ago

This is partially solved. I've yet to figure out what the best mapping strategy is, or how to get the original back. That said ... Someone else dealt with this problem in a different language here: https://stackoverflow.com/questions/51675940/converting-an-image-from-cartesian-to-polar-limb-darkening And this gentleman had a routine to do coordinate mapping in the appendix of his thesis, which conveniently used CImg: https://andrew.ambrose.thurman.org.uk/theses/appendix.pdf

Originally, I wanted to do this using Rcpp, but for whatever reason the example given in the imageR docs, didn't work for me. Specifically, I got a compilation error saying it couldn't find XFreeColormap or something like that. Googling that problem showed that it could be solved by specifying some libraries at the end of the compilation line. But how to do that?

I fixed the that problem by using a plugin, thusly: `library(imager) library(Rcpp)

https://stackoverflow.com/questions/57778952/correctly-registering-a-plugin-to-use-eigen-via-rcpp

Some C++ code as a character string

The function takes an image, applies blur and erode operations, and returns the result

foo.inline <- " NumericVector foo(NumericVector inp) { CImg img = as<CImg >(inp); img.blur(3).erode(4); return wrap(img); } "

Rcpp::registerPlugin( name = "myplugin", plugin = function(x) { plug <- Rcpp::Rcpp.plugin.maker(libs="-lm -lpthread -lX11") settings <- plug() settings } )

Rcpp::cppFunction(code = foo.inline, depends = "imager", plugins = "myplugin", verbose=TRUE,rebuild=TRUE,showOutput=TRUE) That made "foo" available. But I wanted to keep what I was going to be doing in a separate file. So I created a .cpp file to hold this: #include

include "/home/matt/R/x86_64-pc-linux-gnu-library/3.6/imager/include/imager.h"

include "/home/matt/Downloads/imager-master/src/wrappers_cimglist.h"

include "/home/matt/CImg/CImg.h"

using namespace Rcpp; using namespace cimg_library;

// [[Rcpp::plugins(myplugin)]] // [[Rcpp::export]] NumericVector footoo(NumericVector inp) { float theta; CImg img = as<CImg >(inp); theta = img.height(); img.blur(3); img.erode(4); return wrap(img); } Some of the includes probably aren't necessary, but what the heck. By then calling sourceCpp("mattfns.cpp",rebuild=TRUE,cleanupCacheDir=TRUE,showOutput=TRUE,verbose=TRUE) ` "footoo" is now available.

Okay, so now I can compile C++ code and access the CImg API. Moving on ...

I was able to reproduce (pretty closely, I think) the first example on the above linked page, which showed output from openCV. The "polar" routine I wrote was adapted from the thesis noted above (hat-tip to that guy, as well as to the source code examples in the imager-master repository, which showed how to pass and access the image data):

`// [[Rcpp::plugins(myplugin)]] // [[Rcpp::export]] NumericVector polar(NumericVector im, double mfactor=0) { CId imgin = as(im); int imheight = imgin.height(); int imwidth = imgin.width(); int imdepth = imgin.depth(); int imspect = imgin.spectrum();

float theta,radius; float multfactor; int polarx,polary,rectx,recty,expdradius;

if (mfactor==0){ //like stackoverflow/opencv result multfactor = (imwidth/sqrt(imwidthimwidth+imheightimheight)); } else { //multfactor = 1.0; //1=accurate pixels from center? //multfactor = 1.0/2.0f; //=percent from center accurate multfactor = mfactor; } CId imgout(imwidth,imheight,imdepth,imspect); //expdradius = ceil(sqrt((imwidth^2)+(imheight^2))); //CId imgout(expdradius,expdradius,imdepth,imspect); imgout.fill(0);

for (polary=0; polary<imheight; polary++){ theta=(polary1.0f/imwidth)2.0fcimg::PI; for (polarx=0; polarx<imwidth; polarx++){ radius=polarxmultfactor; rectx=(int)(radiuscos(theta)+imwidth/2.0f); recty=(int)(radiussin(theta)+imheight/2.0f); if(rectx<imheight && recty<imwidth && rectx>0 && recty>0){ imgout(polarx,polary,0)=imgin(rectx,recty,0); imgout(polarx,polary,1)=imgin(rectx,recty,1); imgout(polarx,polary,2)=imgin(rectx,recty,2); } } } return wrap(imgout); } The second example in the link above looks different though, and comes from ImageMagick. The convert/distort/polar/depolar routines in IM don't obviously seem available to me in R (would appreciate hearing if someone knows differently) and it that example has the advantage of fulling wrapping the image. So I wanted to reproduce that. The code I came up with is // [[Rcpp::plugins(myplugin)]] // [[Rcpp::export]] NumericVector polar3(NumericVector im, double ctrx=0, double ctry=0) { CId imgin = as(im); int imheight = imgin.height(); int imwidth = imgin.width(); int imdepth = imgin.depth(); int imspect = imgin.spectrum();

int expdwidth = ceil((imwidth/2)*cimg::PI); int expdheight = ceil(imheight/2); CId imgout(expdwidth,expdheight,imdepth,imspect);

int rectx = 0, recty = 0; double indx = 0, indy = 0; float theta = 0.0f;

cimg_forXY(imgout,x,y){ indx = ((x+ctrx)>expdwidth) ? ((x+ctrx)-expdwidth) : (x+ctrx); indy = ((y+ctry)>expdheight) ? ((y+ctry)-expdheight) : (y+ctry); theta = indx2.0fcimg::PI/expdwidth; rectx=(int)(indycos(theta)+imwidth/2.0f); recty=(int)(indysin(theta)+imheight/2.0f); if (imspect==1) { //assume we're dealing with grayscale or tricolor imgout(x,y) = imgin(rectx,recty); } else { imgout(x,y,0)=imgin(rectx,recty,0); imgout(x,y,1)=imgin(rectx,recty,1); imgout(x,y,2)=imgin(rectx,recty,2);
} } return wrap(imgout); } It's called "polar3" for absolutely no good reason. In order to get a close match to the output from IM, I had to use the following call in R: imnew <- mirror(polar3(im,ctrx=604.2,ctry=0.5),"x")`

I don't know if there's a way to change the direction of iterating in Rcpp, so that the subsequent "mirror" operation is not needed. Nor do I have the foggiest idea why I needed to shift the image slightly in order to get a good alignment with the image in the above link. If anyone has actually read this far and has some advice, I'd welcome it.

As it stands now, I'd like to see if there's a better way of mapping the original data to the projection, with any interpolation taking place there. I think but am not sure that it might get rid of some blurriness. But maybe there's no better way.

And it would be helpful to be able to reverse all this to get the original back. But I'm not sure how realistic that is (ie if this is inherently a one-way transform).

Hope this helps someone, mconsidine