tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

whic will fail to find a solution - but at least fails *quickly*, even when given a "plausible" starting point.

Graphing the functions involved it would seem that Re(zIT)=50 can only be obtained for small, NEGATIVE, vaalues of C: somewhere around -1e-10. However, there is no value of 'r' in the 'LineMatch' equation, which will generate, negative values for 'C'.

See the graphs I have added to your worksheet, attached

parasitic2.mw

This post (and the very related post at

https://www.mapleprimes.com/questions/223134-How-Do-I-Use-Maple-To-Remove-Noise-From

asked the basic question - Can Maple be used to perform filtering of images to remove various artefacts such as high-frequency noise, or periodic distortions (or indeed anything else)

My response to both of these threads was - obviously, yes. The required sequence of operations being

  1. Perform a Discrete Fourier Transform
  2. Appropriately filter the 2D spectral data
  3. Perform an Inverse Discrete Fourier Transform on the filtered 2D spectrum

My problem was that each of these "intuitively obvious" operations, involves choices - one has to get involved with "grubby" details, many of which can be confusing. I attempted to produce a "solution" which "simplified" most of these choices, just to show that the required solution could be obtained.

Could I make the solutoin more "efficient"- yes

Would the resuting "efficient" solutoin be more difficult for the OP to follow - yes

So when I answer questions like this, should I use a solutiion which is reasonably efficient, and which the OP has a chance of understanding.: or maybe I shoudl get really efficient (aka non-understandable); generally I will always go for the former

depending on whether you are seeking an "analytic" or a "numerical" solution.

Since you don't post the PDE system, then we have no way of guessing

When you have multiple unit definitions for the same dimension - as in lengths may be expressed in inches, feet, yards, roods, furlongs, whatever - whihc unit do you expect Maple to use for the result of a calculation? - and why would it pick that unit.

Luckily you are given a convert() function: this means that you can convert 10*Unit('feet')+6*Unit('inch') to an expression in units of 'feet' or one  in units of 'inch'. (Or units of 'roods', if you add these to the unit system I supplied earlier).

Your second complaint is essentially the same as the first - why would Maple convert inch^2 to feet^2 unless you asked it to?

Both of these problems are addressed in the attached

unitProb2.mw

if you don't to be bothere figuring out plot domains - then just keep winding up numpoints option. For you latest example, this look pretty good and only takes a second or so

restart;
f := unapply(piecewise(0 <= x and x <= 2 and x-2 <= y and y <= x, y^2+x, undefined), x, y);
plot3d(f, 0 .. 2, -2 .. 2, orientation = [220, 75], numpoints=100000);

 

The big green up-arrow in the MaplePrimes to0llbar is your friend, because there is no prospect that I (or anyone else here) is going to retype your code from a picture, to figure out what is going wrong

As a general rule I admire the courage and persistence of those who freely choose to use Maple's 2D-input mode.

The following works perfectly in 1D input mode


 

"(a)(b";

"(a)(b"

(1)

Copy the string, paste into a new cell, delete and retype the opening quote.

"(a)(`b";

 

"(a)(`b"

(2)

Copy the string from the first cell, add f() around it.

f("(a)(`b");

 

f("(a)(`b")

(3)

 


 

Download 2dIP_again.mw

provided you actually understand the requirement - and I appreciate that this can be tricky

A caveat: I am no kind of expert in 2D image filtering - a very long time ago, I did undersand the principles, but recalling these from still-functioning parts of my wetware is causing brain-ache.

############################################################################

None of the stages I outlined in my explanation is particularly "difficult" to code.

The "difficult" part is determining the nature of the "specific" filter which will remove the undesired artefacts in any given image, whilst causing least possible distortion of the wanted data. Such an "optimal" filter depends not only on the nature of the artefacts which one wishes to remove, but also on the data in a specifc image which one wishes to retain

If the unwanted "lines" in the original image were at 45 degrees, then the "bright" spots in the 2D spectral representation would be  radially symmetric with respect to the centre of the (circularly-shifted) 2D spectrum. A circular band-stop filter would remove these "brightspots". Since the undesired lines are at some random angle, one probably need an elliptical band-stop in the spectral representation, which goes through all four off-centre "bright-spots". One would probably start by using a "brick-wall" ellipse, and having obtained this, "soften" the ellipse so that its cross-section is more "gaussian" than "brick-wall"

I admit that this is finicky to program and is best approached in simple steps (aka don't run before you can walk). So forget everything about image processing and conside how you would achieve the following

  1. Populate a biggish matrix full of ones
  2. For some radius parameter, construct a "circle"of zeros in this matrix. You now have an absolutely brutal band-stop filter which will take out any periodic (45 degree) data.in an image
  3. Now give yourself three parameters and construct an "ellipse" of zeros in the matrix - you can think of the parameters as major axis length, minor axis length, and angle determining orientation of major axis. This too is a pretty brutal band-stop, but with the correct choice of parameters, yu now ought to be able eliminate periodic dat at any angle
  4. Extend both (2) and (3) above to generate an "annulus" of zeros, rather than a single ring. How wide should this annulus be? Another tunable parameter, I'm afraid
  5. "Soften" the annuii generated in (4) above so that instead of a 110000011-type cross-section, it "smoothly"  transitions from 1 to 0 and back again

If you can achieve the above, and more importantly, understand why, you are well on the way to beoming an image filtering guru

@Christopher2222 

You have to be very careful about the distinction between "masking" and "filtering". Consider the steps in the process.

  1. We start with a 2D "image" which I'll call f(x,y) - Technically, ths is a "sampled" image, because 'x' and 'y' are not continuous variables, but represented by discrete entries in a matrix.
  2. Perform a 2D DFT, to created F(u,v) which is a (sampled - see above) 2D spectrum. This will be a matrix of the same size as the input, but will be full of difficult-to-interpret, because it is full of comp[lex numbers
  3. In order to "visualise" this data, I (like most other people) will do the following transformations on F(u,v)
    1. When Maple computes F(u,v), low (spatial) frequencies are in the are represented be entries towards the "edges", and in the "corners" of the F(u,v)-array. High (spatial) frequencies are repesented by entries towards the "centre" of the array. Use of the CircularShift() function means "low" (spatial) frequencies will now be represented be entries along the "central" lines of the array and "high" (spatial) frequencies will be represented by entries towards the "edges" and "corners" of the array. Call this representation CS(F(u,v)), just so that you remeber that it is not the 2-D DFT originally generated
    2. CS(F(u,v)) is still an array of difficult-to-interpret complex numbers, so the next step is compute the magnitude of all the entries, as in abs~( CS(F(u,v)) ). Maybe follow this by some "scaling" of these magnitudes to fit these to a suitable greyscale. The moment one has performed the abs() function, one has ensured that the original data, F(u,v), cannot be recovered from abs~( CS(F(u,v)) ). Furthermore, there will be no way, to recover the original spatial data f(x,y) from abs~( CS(F(u,v)) ).
    3. The only reason for performing these visualisation steps is that it makes it simpler (for mere mortals) to identify regular patterns or anomalies in the data, which one might want to adjust/remove.
  4. One can  now play around with "masks" to see how best to perform the necessary adjustment/removal of data in abs~( CS(F(u,v)) ). to get a mask array M(), with which one can elementwise multiply as in M *~  abs~( CS(F(u,v)) ), to get the "masked" representation. Once one has developed an appropriate "mask", then instead of applying it to abs~( CS(F(u,v)) ), one can apply it to the original 2D spectrum. At this stage one has to circularly "unshift" the mask function to get  CUS(M) which can then be elementwise applied to the original spectral data, as in CUS(M)*~( F(u,v) ). This operation is now technically not masking, but 2D spectral filtering.
  5. Consider a couple of "simple" examples (refer to 3.1 above, for locations of high/low spatial frequencies in the orginal data)
    1. "Brick wall" low pass filtering. Since high (spatial) frequencies are concentrated in the central entries of F(u,v). any "mask" array consisting of an "annulus" of ones around the edges, and a block of zeroes in the central area, when multiplied elementwise with the original data will perform a brick-wall low-pass filter. It is informative to apply the CircularShift() trick to this "mask" to see what it would "look" like at stage (4) above. This will have a central, vertical band of ones, a central, horizontal band of ones, and blocks of zeroes in all four corners
    2. "Brick wall" high-pass filtering. By the same argument as above the "mask" array will have a block of ones in the central area, and with an annulus of zeroes. Again doing the CircularShift() on this array is informative when considering stage (4) above
    3. Consideration of the above two examples, should convince you that a "brick wall" bandpass wilter would be given by some sort of (probably) circular "annulus" of ones, with a block of zeros in the centre, and a further "annulus" of zeroes around the edges.
    4. Obviously one would rarely use this kind of "brick wal" filtering on the 2D DFT data, because as well as filtering, it may/will introduce some odd artefacts when the the original image is reconstructed. It would be more common to use something like a smoothly-varying 2D gaussian.expression. With this giving values of one in the centre of the mask array, tapering to (near) zero along the edges, you have a "smooth" high-pass filter. If you flip this gaussian so that it is giving (near) zero in the corners and one in the centre, you have a "smooth" high-pass filter. Similarly, producing a gaussian-weighted annulus will give a smooth bandpass

Determining the "best" filter to perform a specific operation on a given image is a bit of a "black art", although I am assured that if you work in the area of image processing, then you get quite good at making "educated guesses", with or without the visualisation stages at (4) above

 

You state that you couldn't get " the command from DiscreteTransforms to work". This is worryiing and may be a Maple version issue. I have checkd the following in Maple 18, Maple 2015, Maple2016 and Maple 2017. It works in all of them. If it doesn't work for you please specify the Maple version you are running


 

#
# Initialise.
# Load the ImageTools package
# Load the DiscreteTransforms package
# Load the ArrayTools Package
#
  restart;
  with(ImageTools):
  with(DiscreteTransforms):
  with(ArrayTools):
#
# Read and view the loaded image. Note that the file path/name
# for the image source is only relevant on my machine. OP may
# have to adjust the either/both
#
  img := Read("C:/Users/TomLeslie/DeskTop/clown-noise.jpg"):
  Embed( Scale(img,4)):
#
# Compute the 2D Fourier transform of the loaded image. I'm
# sure that "Maxim" recommended the use of the SignalProcessing
# package commands DFT of FFT. However my reading of the help
# for the SignalProcesing commands is that they are only usable
# for 1D signals. Rather than have to debug this possibility
# I recommended using the DiscreteTransforms:-FourierTransform()
# command. This will use the FFT algorithm whenever the pixel
# count is a power of 2: it will use an efficient DFT algorithm
# whenver the pixel count is composite number (ideally with a small
# number of small prime factors); and if neither of these apply
# then it will perform a DFT anyway; although when one has to do
# the last of these then execution times may start to climb
#
  img2:=FourierTransform(img):

#
# Note that the output of the FourierTransform() command is a
# 128x128 complex array. If one wants to "visualise" this data
# then the conventional approach woul be to compute the absolute
# value of each pixel and then the "log" of the resulting value,
# cos our eys, (like our ears) have a logarithmic response.
#
  img2_vis:= FitIntensity
             ( Create
               ( ln~(abs~(img2))
               )
             ):
  Embed
  ( Scale
    ( CircularShift
      ( img2_vis,
        Height(img2_vis)/2,
        Width(img2_vis)/2
      ),
      4
    )
  );

 

 


 

Download imageProc.mw

The attached i a variation on the worksheet I posted earlier. I have added an execution group which uses the value method returned by pdsolve(..numeric) and used this to compute numerical values for all parameters fo a variety of arguments - if thhis isn't what you want need, let me know

pdeSol4.mw

then

`*`(1/x,y,z);

will work.

I have no idea what you are hoping to achieve with the 'map' function. Form the first lien in the help page for 'map' it says (my emphasis)

apply a procedure to each operand of an expression

Now how, exactly would you expect a binary operator (ie requires two arguments), such as '*' to be applied to " each operand of an expression"

This is essentially a 2-D signal processing problem, which in turn can be reduced to a matrix manipulation problem, so in principle Maple can do it.

From the pictures/spectra which you supply, the periodic "noise", ie lines which appear in the spatial domain, will (because they are strictly periodic) be converted to points in the frequency domain. If the lines were perfectly vertical, then this would show as couple of points on the horzontal axis in the frequency domain. If the distorting lines were perfectly horizontal, this would show as a couple of points on the vertical axis in the spectrum. Since, in your example, the lines are at an angle, then you will get a mixture of the two limiting cases, which would be four points, symmetriccwrt the origin, but not located on either axis..

Since these lines are at a very specific frequency, then I would have thought that the optimal way to remove them would be with some kind of 2D bandpass filter

Failing the above, I would probably start with

#
# Generate some more or less random plot
#
   p := plots:-polyhedraplot( [0,0,0],
                                                polytype=SquareOrthobicupola
                                             );
#
# Export this plot to s suitable location. Note that OP will
# have to use some suitable file path for his/her machine
# rather than
#
#  "C:/Users/TomLeslie/Desktop/
#
    Export("C:/Users/TomLeslie/Desktop/figure.wmf", p);

From within MS Word, use Insert->Picture, and access the file name used in the above "Export" command. If this produce some loss of  quality, please provide an example plot so I/we can check

Multiple problems

  1. There is no package "OrthogonalExpansions" in Maple 2015, so I assume that you have a Maple "add-on" installed to which I have no access
  2. There is no command "FourierSeries" in Maple 2015, so I assume that this command is defined within the "OrthogonalSeries" add-on package above.
  3. GIven that I have no idea what add-on package you have loaded, and therefore cannot execute your worksheet, my ability to solve your problem is limited.

The general method for  "animating" a list of plots simultaneously, is something like (warning - "toy" example)

P := animate(plot, [ [seq(sin(k*x+t),k=1..4) ],x=-Pi..Pi ], t=-Pi..Pi, frames=24):
plots:-display(P, insequence=true);

where [ [seq(sin(k*x+t),k=1..4) ],x=-Pi..Pi ], defines the list of plots you want to "animate", t=-Pi..Pi, defines the animation variable, and frames=24 defines the number of "frames" in the animation

 

 

 

First 109 110 111 112 113 114 115 Last Page 111 of 207